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[57] ABSTRACT 

An ophthalmic diagnostic instrument determines the 
shape of the cornea through projection of an image onto 
the cornea through the optics of the diagnostic instru- 
ment. The instrument and the method of the invention 
involve folding a projected pattern of discrete separated 
point light sources so that the pattern is projected 
toward the eye coaxially with return collected light 
reflected off the cornea. The instrument avoids any 
need for a pattern light source directly adjacent to the 
eye, and provides the surgeon or other eye care special- 
ist with a real time image accurately displaying the 
shape of the cornea. The surgeon is thus able to monitor 
the corneal shape prior to surgery, to monitor its 
changes during the course of the surgery, and to further 
monitor the cornea in post operative stages. In a specific 
embodiment of the invention, a real image of the pattern 
of point light sources is formed inside or very closely in 
front of the objective lens of the system so that the 
objective lens becomes a field lens and the angle of view 
of the system is enlarged. 

30 Claims, 8 Drawing Sheets 
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OPHTHALMIC DIAGNOSTIC APPARATUS AND 
METHOD 

BACKGROUND OF THE INVENTION 

This invention relates to ophthalmic analytical and 
diagnostic systems, and in particular the invention is 
concerned with obtaining accurate determinations of 
the shape of the human eye structures such as the cor- 
nea, the lens and the retina. An apparatus in accordance 
with the invention measures, calculates and displays the 
shape of selected cross sections of the cornea, for exam- 
ple, and is intended for use by ophthalmic surgeons as 
well as the eye care community at large- 
One of the principal activities of the eye care special- 15 
ist. which includes both ophthalmologists and optome- 
trists, is to determine the refractive power of the eye as 
an optical system. Since the only major refractive index 
change along a light path entering the eye to impinge on 
the retina occurs on the first air to tear layer interface 20 
or, approximately, at the corneal anterior epithelial 
surface, the precise measurement of the shape of the 
corneal epithelial surface is the key to estimating the 
refractive power of a given eye. 

Traditionally, the eye care specialist has been satisfied - 5 
with a measurement from keratometric readings. The 
keratometric readings ( 4 * K-readings") correspond to the 
curvature of the corneal epithelial surface at the inter- 
section of the corneal epithelial surface with the central 
visual axis of the eye. The K-readings are usually dis- 30 
played in diopter power which is proportional to the 
reciprocal of the radius of curvature. The K-readings 
provided by keratometers correspond to the curvatures 
at one point on the corneal epithelial surfaces along two 
surface rays passing through that point. Usually, the 35 
two rays match the semi-major and semi-minor axes of 
the eye which are the nasal-temporal (horizontal) axis 
and the superior-inferior (vertical) axis of the eye. 

Since the first concern of the eye care specialist is 
central, axial vision, the K-readings, which only pro- 40 
vide two curvature measurements along the semi-major 
and semi-minor axis normal to the visual axis, represent 
a fair estimate of refractive power along the most criti- 
cal light paths in the human eye. 

To date, the eye care community has relied on the eye 45 
surface being approximated as a combination of a 
sphere and a cylinder, thus the reference to 20/20 as a 
visual standard. This approximation is exact at the inter- 
section of the corneal anterior surface with the visual 
axis. The approximation is known to fail as one pro- 50 
ceeds radially outward from the central visual axis 
towards the limbus, roughly the outermost edge of the 
cornea where the triple-point transitions between cor- 
nea, sclera, and iris tissues take place. 

There are several instruments for measuring the loca- 55 
tion of the corneal anterior surface in proximity of the 
limbus as well as in the central region, but the display 
generated from these measurements usually assumes 
that the eye can also be approximated as a combination 
of spheres and cylinders. In a sense, these instruments 60 
spread the error in approximating the shape of the cor- 
neal surface from being concentrated toward the limbus 
to being distributed over a greater region. 

Notable exceptions are instruments based on confocal 
microscopy that measure the actual curvatures without 65 
simplifying assumptions. However, systems based on 
confocal microscopy have very limited fields of view, 
considerably smaller than the full corneal surface. Such 



systems must then rely on a sequence of measurements 
over time which are subsequently made to piece to- 
gether using either fractal techniques or some boundary 
matching algorithm. These paste-ups involve some 
form of interpolation, albeit on the boundary of the 
images rather than in the interior. In contrast with inter- 
polations which are based on a sparse set of measure- 
ments, confocal techniques provide dense measure- 
ments at the expense of not having them performed 
simultaneously. Even though the sequential measure- 
ments can be formed very quickly, involuntary eye 
motion is known to occur within millisecond time 
scales, faster than the time required to complete data 
gathering using confocal techniques. This can introduce 
errors. 

Since only a finite number of measurements of the 
actual location of the corneal anterior surface are possi- 
ble, interpolation techniques are an intrinsic part of 
displaying a continuous shape based on the measured 
information. In most instruments still in use by the eye 
care community, the error in measuring the shape of the 
corneal anterior surface is often not in the measurement 
technique but in the numerical interpolation techniques 
utilized to prepare the display of the continuous corneal 
surface cross-section. 

SUMMARY OF THE INVENTION 

In accordance with the present invention, a system, 
apparatus and method for corneal shape determination 
and other eye structures achieves a high degree of accu- 
racy by real-time measurements, calculations and dis- 
play. 

An interpolation technique used in accordance with 
the invention is based upon not only measuring the 
location of a given set of points on the surface of the 
corneal tear layer, but in performing the measurement 
in such a manner as to be able to solve the non-linear 
ordinary differential equation describing the surface in 
real time. This technique further provides several 
higher order derivatives which are used in generating 
the continuous corneal shape. 

Aside from improving the K-readings while simulta- 
neously reducing the global error in approximating the 
shape of the cornea's anterior surface, the process and 
apparatus of this invention are designed to provide the 
actual corneal shape in real time. By real time is meant 
faster than an eye care practitioner's capacity to observe 
and recognize a change in corneal shape, generally on 
the order of one second or less. 

Several instruments have claimed capability of mea- 
suring a corneal cross-section shape in real time, as 
defined herein. These have relied on a combination of 
scanning techniques, confocal microscopy techniques, 
and image reconstruction techniques. Namely, using 
confocal microscopy, a high-resolution set of measure- 
ments is taken of selected small segments of the corneal 
cross section and the shape of the individual crosssec- 
tion pieces is analyzed and then fitted together with the 
other separately scanned shapes to provide an approxi- 
mating shape. 

Given that the human eye is constantly in motion 
either through voluntary or involuntary actions and 
that these motions do not necessarily correspond to 
rigid body motions, scanning techniques which later 
patch up different pieces of the corneal cross section 
may have errors introduced by eye motions and eye 
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distortions. The errors may be diminished by utilizing 
high-speed scanners, but are not removed. 

The advantage of the described scanning technique is 
that it allows for performing a high number of measure- 
ments in a small region, thus generating very high reso- 5 
lution and accurate definition of the shapes in the lim- 
ited area observed. To achieve the same resolution 
globally, that is from limbus to limbus, would require an 
unwieldy instrument with many light sources and mea- 
surement points. Hence speed is achieved with the scan- 10 
ning technique as well as high local resolution at the 
expense of global uncertainty due to the patchwork 
effort. 

In one of the embodiments of the present invention, 
the instrument achieves comparable high local resolu- 15 
tion without risking the uncertainty of global accuracy. 
This is achieved by performing a global measurement 
whereby the field of view of the instrument is adjust- 
able. Thus, when global information is desired as to the 
shape of the cornea, one measurement rather than a 20 
sequence of scans is performed. And when a high reso- 
lution measurement is required of a particular region of 
the corneal surface, one of the embodiments of the 
invention will be able to narrow the same number of 
measurements into a limited field" surrounding the de- 25 
sired corneal region. 

The approach is comparable to the use of photoli- 
thography to generate pattern cuts in the semiconduc- 
tor industry in the sense that a large template is progres- 
sively focused into smaller regions to generate tighter 30 
and smaller effects. In preferred embodiments of the 
present invention a zooming technique is used to define 
the field of view, and a corresponding zooming tech- 
nique is used to enlarge the field to fill the display moni- 
tor. 35 

One of the important considerations to ophthalmolo- 
gists seeking to perform corneal reconstructive surgery, 
whether by radial keratotomy, corneal epikeratoplasty, 
keratimileusis, keratophakia, wide area laser ablation, or 
other procedures, is to accurately and reliably measure 40 
the shape of the cornea. This is important not only prior 
to the initiation of a surgical procedure, but also during 
evolution of the shape in the surgical procedure, and 
after the surgery, as the healing process takes over post- 
operatively. 45 

One of the features of the present invention is not 
only to satisfy the corneal global measurement needs of 
pre- and postoperative surgical procedures, but to like- 
wise provide the high resolution needed to follow the 
effects of surgery during the course of the procedure. 50 
For example, to determine the depth of an incision 
during a radial keratotomy procedure, the surgeon pre- 
selects the depth of the protruding diamond blade from 
the scalpel and depresses the full depth of the blade into 
the cornea. Surgeons are therefore relying that a uni- 55 
form thrust pressure on the blade point can be main- 
tained as the scalpel is traversed over the cornea and 
that the shape of the cornea does not deform during the 
course of the procedure. It has been often observed that 
the cornea deforms during the course of radial keratot- 60 
omy. Thus, the surgeon is left guessing and relying on 
his own intuition as to the incision depth of radial cuts. 
The depth of the incisions becomes progressively un- 
certain with each successive cut because of the in- 
creased deformation. 65 

In one of the embodiments of the invention, the sur- 
geon is enabled to first determine the corneal shape 
. immediately prior to commencing the procedure, then 



to observe the measured incision depth of each cut as it 
is being performed, to readjust his blade progressively 
as needed based upon actual incision depth measure- 
ment, and then once again to provide a global measure- 
ment of the resulting corneal shape. 

Another problem addressed by the present invention 
again involves the utility of a measurement device dur- 
ing surgery. Keratometers have currently been used in 
surgical theaters, such as the Terry keratometer. but 
their efficacy has been limited not only by the extent of 
the information provided, but by their accuracy and by 
the obstruction of the patient's eye to the surgeon while 
the measurement is being performed. One of the reasons 
for this obstruction is mandated by the need to provide 
a field of illuminated points which the keratometer then 
detects as reflections from the cornea. 

In order to get reflections or data points near the 
central visual axis, keratometers have needed to place 
such illumination points near the central visual axis. . 
These illuminators are bulky and get in the way of the 
surgeon's access to the eye on which he is operating. In 
one of the embodiments of the present invention, this 
problem is solved by placing the illuminators a consid- 
erable distance away from the eye, folding an image of 
the illuminators into either a surgical microscope or 
other imaging apparatus via a beam splitter, and then 
rather than physically placing the illuminators along the 
visual axis, the present invention projects a real image 
of the illumination points at the location where the 
illuminators would have been required. This location is 
between the patient's eye and an objective lens of the 
instrument, or of a surgical microscope to which the 
instrument is attached. The real image is reflected off 
the surface or surfaces of interest in the eye, and re- 
flected illumination points parallel to the optical axis of 
the instrument are collected and detected through the 
instrument. A real-time display is generated, preferably 
with ocular cross sections as selected by the surgeon, 
along with numerical topological data. 

An important aspect of the present invention is that 
the optics of the system use the objective lens as a field 
lens for the pattern image and that the optics relay the 
Fourier plane of the objective lens, located behind the 
objective lens in the system, to a relayed, distant posi- 
tion in the instrument. This gives the opportunity and 
the spatial distance to fold in one or two light source 
patterns, between the Fourier plane of the objective 
lens and the relayed or transferred position of the Fou- 
rier plane. In one embodiment, as discussed below, two 
different light source patterns may be folded into the 
optical axis of the instrument, using two different beam 
splitters. 

In this regard, two separate displays can be formed 
for real time review by the surgeon. One, a qualitative 
image showing elevation contours of the eye, derived 
from a series of concentric light rings; and the other 
from a selected pattern of discrete points of light, for 
quantitative analysis in producing cross sectional repre- 
sentations of the shape of the cornea. 

The utility of the invention is not restricted to im- 
proving radial keratotomy procedures. Any surgical 
procedure which seeks to alter the refractive power of 
the eye benefits from having accurate displays showing 
the course and effect of the procedure. More generally, 
any surgical procedure which invades the eye and 
which in turn necessitates wound closure can be greatly 
benefited by the process and apparatus of the invention. 
Also, the instrument can be used for purely diagnostic 
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purposes, such as by an optometrist for fitting contact 
lenses. 

The various embodiments described correspond to 
different configurations depending on the actual needs 
of the surgeon. The common feature is to provide high 
resolution wherever it is mandated while preserving 
computational speed and global accuracy wherever 
extreme resolution can be relaxed. 

It is therefore among the objects of the invention to 
enhance ophthalmic diagnostic and surgical procedures 
by providing an apparatus, system and method for high 
speed, real time precision monitoring of the shape of the 
cornea, both epithelial and endothelial surfaces, and of 
other ocular surfaces. These and other objects, advan- 
tages and features of the invention will be apparent 
from the following description of preferred embodi- 
ments, considered along with the accompanying draw- 
ings. 

DESCRIPTION OF THE DRAWINGS 

FIG. 1 is a schematic diagram showing the layout of 
an ocular diagnostic system and apparatus in accor- 
dance with a preferred embodiment of the present in- 
vention. The figure shows a series of optical elements 
which can be incorporated in this embodiment of the 
invention. 

FIG. 1A is a schematic diagram showing in greater 
detail certain portions of the system of FIG. 1. 

FIG. IB is a view similar to FIG. 1, but showing a 
slightly modified form of the system, and with an exam- 
ple of distances and other optical values which can be 
used in the system. 

FIG. 2 is a schematic diagram showing portions of 
the apparatus of FIG. 1 as they can be incorporated into 
a surgical microscope, preferably with a simple auxil- 
iary camera mount connection (for example, a C-mount 
connection). 

FIG. 3 .is a schematic view showing an example of a 
video display which can be presented as a result of the 
information gathered by the instrument of the inven- 
tion. 

FIG. 4 is a view showing an example of quantitative 
information which can be displayed to the user of the 
instrument. 
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FIG. 11 is a schematic view showing a reflected 
pattern of discrete point light sources which might be 
received from the measurement projection shown in 
FIG. 6, and also showing an example of a secondary 
reflection which is taken from the endothelial surface, 
i.e. the back surface of the cornea. 

FIGS. 12. 13 and 14 are schematic graph plottings 
showing examples of light intensity versus distance 
across the cornea, for a reflected projection such as 
shown in FIG. 11. and indicating analysis of these re- 
flections to obtain information as to both the front sur- 
face and the back surface of the cornea, i.e. the epithe- 
lial and endothelial surfaces. 

DESCRIPTION OF THE PREFERRED 
EMBODIMENTS 

In the drawings. FIG. 1 shows in schematic represen- 
tation a system of optical elements in accordance with 
the invention for use in carrying out ophthalmic diagno- 

20 sis and analysis. 

The system, generally identified by the reference 
number 10, includes an illuminator or light source 12. a 
pattern plate or disk 14 having a pattern of holes cut in 
the plate for producing a desired pattern of discrete 

25 light sources, a non-distorting beam splitter 16, a lens 18 
which projects an image of target 14 onto an image 
plane at 22. This image plane 22 is close to or coinciden- 
tal with the system of objective lens 20. The purpose of 
placing the image at this location 22 is to have the ob- 

30 jective lens 20 serve as a field lens, that is bending the 
rays of light that form the image towards the patient's 
cornea 24. 

As indicated in FIG. 1, the focused image 22 of the 
pattern is a real image, formed at some plane at or near 

35 the lens 20 and between the lens and the patient. The 
real image preferably is in the lens 20, but it can be very 
closely in front of the lens (i.e. a few millimeters in 
front). In this real image, each point source of light 22a 
projects a cone of light toward the patient. Thus, each 

40 point source 22a in the real image makes an infinite 
number of specular reflections off the front surface of 
the cornea 24 of a patient's eye 26. As explained below, 
the F-number of the final lens 20 determines the maxi- 
mum area of the cornea that can be measured. The 



FIG. 5 is a view showing a conventional surgical 45 objective lens serves as a field lens, and the patient's 



microscope with which the embodiment of FIG. 2 can 
be used. 

FIG. 6 is a schematic view showing a simple pattern 
of rectilinear sequences of point light sources which can 
be used for measuring the cornea in accordance with 50 
one embodiment of the invention. 

FIG. 7 is a schematic view similar to FIG. 6, showing 
another pattern which can be used for measuring the 
cornea. 

FIG. 8 is another schematic view showing concentric 55 
circles of light which can be projected on the cornea 
simultaneously with the pattern shown in FIG. 6 or 
FIG. 7, strictly for qualitative information for the sur- 
geon, for comparing against the measurements deter- 
mined from the quantitative measurements obtained via 60 
the pattern of FIG. 6 or FIG. 7. 

FIG. 9 is a schematic view showing a distortion of the 
pattern shown in FIG. 7, as an example of what may be 
read and analyzed by the apparatus of the invention for 
determining cornea shape. 65 

FIG. 10 is another schematic view, showing a re- 
flected pattern as produced by the projection shown in 
FIG. 8, and the distortion of the reflected pattern. 



cornea must be at the focal length of the lens 20. This 
assures that the light reflected off the eye parallel to the 
optical axis of the instrument is then brought to a point 
behind the lens 20 at the focal distance of the lens 20. 
This enables the return light to be apertured down as 
discussed below, to select only those rays which were 
paraxial off the eye. This enables the system to localize 
a detected point to a point on the cornea from which 
that ray was reflected. If the objective lens 20 were not 
situated to serve as a field lens, outermost points of light 
in the pattern would not reflect off the cornea. As a field 
lens, the lens 20 efficiently bends the outer points of 
light toward the eye. 

It is preferred that the focal length of the lens be great 
enough to provide an unobstructed, comfortable dis- 
tance from the instrument to the patient and adequate 
working room for the surgeon, for surgical applications. 

The F-number of the objective lens 20 is most impor- 
tant in its function as a field lens as it will determine the 
maximum angle from the optical axis at which a ray can 
be reflected from the cornea parallel to the axis. If for 
example a commercially available F/2 lens is used, then 
the region of coverage will be about 3 mm diameter on 
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the cornea. A lens with a smaller F-number will cover 
a proportionally larger region on the cornea. 

As indicated in FIG. 1A. each of the real-image point 
light sources 22a makes at least one reflection 22b 
which will be parallel to the central axis of the objective 5 
lens 20. with all axial reflected rays 22b parallel as 
shown in the drawing. For each point light source 22a, 
the reflected axial ray 22b will be unique unless the 
corneal surface has extremely strong local imperfec- 
tions or distortions in the corneal curvature, which 10 
could theoretically cause more than one reflected axial 
ray 226 to occur, from spaced locations on the cornea. 

Other rays of light reflected off the cornea will reach 
and pass through the lens 20. but as will be seen below, 
only those returning reflected rays which are very 15 
nearly parallel are passed through the system for analy- 
sis. Those are the rays and points which will supply data 
points to be compared with the original pattern as pro- 
jected through the plate 14 to supply data which can be 
solved to determine the shape of the cornea. 20 

As shown in the overall schematic view of FIG. 1, 
the returning reflected rays pass back through the lens 
20. then through the lens 18 and the beam splitter 16, an 
aperture or spatial filter 30 and a further lens 32, ulti- 
mately to be focused on a detector or camera plane 34. 25 

The curvature of the cornea 24 forms a virtual image 
93 of the target image 22. In the article "Suggested New 
Methods for Photokeratoscopy, a Comparison for 
Their Validities, Part I", by S. G. El Hage, American 
Journal of Optometry and Archives of American Academy 30 
of Optometry, November 1971, El Hage pointed out that 
an aperture or spatial filter at the back focal plane or 
Fourier plane of the objective lens 20 will only pass rays 
parallel to the axis thus localizing those rays from a 
given point of the virtual image 93 to those that are 35 
reflected from a specific point on the cornea 24. In this 
embodiment, it is desired to have space behind the ob- 
jective lens 20, the lens 18 is used to relay the Fourier 
plane of the lens 20 to the aperture 30. The aperture 30, 
being in an image of the Fourier plane, will likewise 40 
select only those rays reflected from the cornea 24 par- 
allel to the axis. 

The rear lens 32 of the system focuses a distorted 
image of the virtual image 93 of point light reflections 
on the detector or camera plane 34. 45 

As shown schematically in FIG. 1A, the camera 
plane 34 has a central axis c which lies on the optical 
axis of the system, including the objective lens. Ideally 
this axis is placed as closely as possible to the center of 
the cornea or visual axis v. If these axes are significantly 50 
displaced, then much of the light reflecting off the cor- 
nea will not be returned through the system. This dis- 
cussion assumes the axes coincide, but adequate infor- 
mation can be obtained over small deviations (e.g. one 
millimeter). If a reflected, returned point lies on the 55 
center axis c of the camera plane 34, then that ray ema- 
nated from the visual axis v of the cornea, at least as 
respects one orthogonal direction on the cornea and on 
the camera plane 34, which is shown as the left-right 
direction in the plane of the paper in FIG. 1A. 60 

Likewise, if a particular point of light is focused on to 
the camera plane or detector face 34 at a distance x' 
from the center axis c, that distance corresponds to, and 
is linearly proportional to, a distance x of the reflecting 
point on the cornea for that ray 22b as measured from 65 
the visual axis v. If a depth distance y is determined, 
measured from an arbitrarily chosen datum d to the 
point of reflection on the cornea, and a series of such x 



and y can be determined, then a differential equation 
can be solved to define y as a function of x. giving the 
curvature of the cornea in this direction or along the 
subject axis, i.e. in the plane of FIG. 1A. Similarly, 
measurements and calculations can be made along an 
orthogonal axis on the cornea (e.g. the nasal-temporal 
and superior-inferior axes can be used), giving as much 
information regarding the cornea's shape as is normally- 
needed for any diagnostic or surgical procedure. 

The y distance indicated in FIG. 1A can be derived 
through information regarding the degree of distortion 
of the reflected point light pattern, and the spatial rela- 
tionship among the points of light, as compared to the 
pattern as originally projected and as arranged in the 
real image 22. Thus, considering the parallel ray 22b in 
FIG. 1A, which is shown as emanating from the real 
image point light source 22a on the right in FIG. 1A. if 
the cornea curvature is less steep at the point of reflec- 
tion, i.e. at a shallower angle with respect to a tangent to 
the cornea at the visual axis, then the parallel ray 22b 
would have originated from a different real image point 
source, one farther to the left in the pattern. The right- 
end point source 22a would have created a parallel 
reflection only from another point on the cornea, far- 
ther to the right as viewed in FIG. 1A. Each of the 
reflected points as detected at the camera plane 34 can 
be identified electronically, essentially by counting 
points in the array. 

FIG. 6 shows one example of a projected light pat- 
tern which can be used in the system and method of the 
invention. In this simple pattern, a vertical rectilinear 
array 40 is crossed orthogonally with a horizontal recti- 
linear array 42, with the intersection point correspond- 
ing to the visual axis of the eye. This is the simple pat- 
tern assumed with reference to FIG. 1A. 

A more complex pattern 44 of points is shown in 
FIG. 7. This pattern, shown as an asterisk-like pattern of 
linear arrays of points, gives data from many more 
points on the cornea. It may define an outline of a five- 
pointed star or any similar type of pattern, but prefera- 
bly it has some means of identifying its rotational orien- 
tation. It may have an outline of a star with an odd 
number of points, so that the asymmetry can help iden- 
tify the detected, reflected points by correlating them 
with the originally projected pattern 44. FIG. 9 shows 
an example of a reflected pattern 46 which might result 
from the pattern 44 shown in FIG. 7, as reflected from 
a cornea with some degree of distortion. 

FIG. 1 schematically indicates that the detector or 
camera plane 34 is connected to a microprocessor 50. 
The microprocessor may be connected to a display 
device, such as a CRT monitor 52 as indicated. Data 
gathered from the system as described is received by the 
microprocessor 50 and analyzed. Each detected point is 
correlated with the location of the particular point in 
the source pattern from which it emanated. The x value 
is determined for each point, i.e. the distance from the 
optical axis v from which the point was reflected off the 
cornea. This is determined by direct proportioning, 
from the known magnification of the system. Each 
reflected point has an x value which is the distance from 
the optical axis of the system. Each linear array of 
points in the image must be separately analyzed and 
fitted to the mathematical approximation. If the com- 
plex pattern 44 shown in FIG. 7 is used, formed of an 
asterisk-like array, the analysis and computation are 
made along each line of the pattern. 
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By the method and system of the invention, the math- 
ematical shape of the cornea is determined by assuming 
an analytical approximation to the surface shape. The 
analytical approximation is then substituted into a dif- 
ferential equation and some type of appropriate fit is 5 
performed to determine the coefficients that satisfy the 
differential equation. In a preferred embodiment of the 
invention, a nonlinear least squares fit is performed. 

These operations are performed in the microproces- 
sor 50. The processor has programming to review a io 
great number of x values as determined on the detector 
34. substituting all of these values into the differential 
equation and arriving at a formula for y and as a func- 
tion of x. 

A differential equation suitable for this purpose is 15 

20 

where y is the depth of the reflection site away from a 
datum plane (such as the datum plane d shown in FIG. 
1A), x is the distance from the visual axis, and a and b 
are coordinates representing the location of the real 
image of the illumination point in space. *a is a distance 25 
of the particular illumination point 22a (see FIG. 1A) 
from the visual axis and b is the depth of that illumina- 
tion point out from the datum plane d. 

The differential equation used in this process is not 
new. It is a general equation which can be used to repre- 30 
sent the shape of any surface, and is described in the 
article "Suggested New Methods for Photokeratos- 
copy, a Comparison for Their Validities. Part F\ by S. 
G. El Hage, American Journal of Optometry and Archives 
of American Academy of Optometry, November 1971, 35 
page 897. In the article, El Hage discusses various uses 
of this general equation for solving the shape of the 
corneal surface. Also, he relates the corneal surface 
shape to one of the keratoscope rings in photokeratos- 
copy. Thus, this derivation in itself does not form a part 40 
of the present invention, but is hereby incorporated in 
this application by reference as illustrating that such 
derivation is known in the art. 

At page 909, El Hage shows an optical arrangement 
for projecting an image onto a cornea and for detecting 45 
reflected light from the cornea. His source is analogous 
to the real image in the present invention, and El Hage 
had a number of optical elements between the source 
and the eye, including a beam splitter between the ob- 
jective lens and the eye. 50 

Returning to FIG. 1, the illuminating light source 12 
may be a visible light source, in preferred embodiments 
of the invention wherein the system is not combined 
with a coaxial surgical laser. For example, an incandes- 
cent lamp can be used. The pattern plate or target 14 55 
may be laser or photolithographically cut, with hole 
sizes on the range of about thirty microns. The beam 
splitter 16 may be a simple nondistorting plate glass 
beam splitter, with a surface coating of about 50% re- 
flectivity. 60 

In one specific embodiment of the invention, particu- 
lar lenses and lens relationships may be selected as indi- 
cated in FIG. IB. In FIG. IB the distances between 
lenses, focal lengths and diameters of the various lenses 
are given for this specific embodiment. Other relation- 65 
ships and distances are also given, including the diame- 
ter of the aperture or spatial filter 30. The system of 
FIG. IB shows a single light source 12 projecting a 
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pattern and being folded into the optical axis of the 
instrument. 

The detector or camera plane may comprise a high 
density photodetector array, for example. 

As indicated schematically in FIG. 1, the micro- 
processor 50 is connected to a display monitor 52. An 
example of the type of display that can be presented to 
the physician in real time is indicated in FIG. 3 by the 
reference number 54. In the upper left quadrant of the 
screen, patient identifying data is given, along with 
K-readings and thickness readings. A more detailed 
example of this information is shown in FIG. 4. 

The lower left and lower right quadrants of the dis- 
play 54 show examples of depth references of the epi- 
thelial and endothelial cornea surfaces at cutting planes 
A and B shown in the plan view of the upper right 
quadrant. The location of these cutting planes is prefer- 
ably selectable by the physician, via inputs to the micro- 
processor 50 (not shown). 

The distorted image 56 shown in the upper right 
quadrant of FIG. 3 is derived from a second projection 
which is preferably included in preferred embodiments 
of the invention. As illustrated in FIG. 1, a second pro- 
jection may be folded onto the axis of the lens system 
via a second beam splitter 60, which reflects light from 
an illuminator light source 60 to a pattern or mask 64. 
The mask 64 has a plurality of concentric circle cuts so 
as to project a real image of the concentric circles in 
front of the cornea as is done with the pattern 22 of 
point light sources. 

FIG. 8 shows schematically a series of concentric 
circles in a pattern 55 which can be projected via the 
pattern plate 64. The detector 34, which may be a pixel 
array of very high density, can receive and detect both 
reflected images simultaneously. The concentric ring 
pattern can be discerned from the point source pattern 
by the contiguity of each ring. The software employed 
by the microprocessor 50 can sample each pixel receiv- 
ing light and determine whether any immediately adja- 
cent pixel is also receiving light. If so, the contiguity of 
a ring is indicated. In contrast, the patterns of point light 
sources such as shown in FIGS. 6 and 7 will not display 
appreciable contiguity. Thus, the microprocessor 50 
can separate these images and analyze each separately. 
Alternatively, in a separate embodiment of the system 
and additional camera detector can be placed together 
with an additional beam splitter to separate the image of 
the continuous rings from the image of the discrete 
point sources. 

As in a conventional corneoscope or in using a Pla- 
cido ring, the concentric light rings produce a reflection 
off the cornea which is distorted in a way correspond- 
ing to distortions on the corneal surface. This can result, 
for example, in a pattern of distortion 56 such as shown 
in FIG. 3. 

FIG. 1 shows that, with two different light patterns 
folded into the system, onto the axis of the lenses 18 and 
20, polarizers 66 and 68 should be used to establish 
opposite polarity for the two different images being 
projected. 

A polarizer used as an analyzer 94 may be rotated to 
select either of the projected images. 

FIG. 2 is a schematic representation of an alternate 
embodiment showing some of the same elements pres- 
ent in the embodiment of FIG. 1, but in an arrangement 
for connection directly with a surgical microscope. 
Surgical microscopes, such as those made by Week, 
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Nikon, Topcon. Zeiss. Nidek. or Wild, usually include a 
standard auxiliary bayonet mount or screw attachment 
for a camera. FIG. 5 shows a typical standard surgical 
microscope. An auxiliary mount 70 (e.g. a C-mount) is 
shown in FIG. 5 and schematically indicated in FIG. 2 5 
as coupling the system embodying the elements 12. 4, 
34, 32, 18 and 16 to a fitting or optical tube 72 on the 
surgical microscope. Generally the surgical microscope 
will have optics to produce an image at an image plane 
74 which is a standard distance from the auxiliary *° 
mount on the fitting 72, for coupling a standard 35 milli- 
meter video camera to the surgical microscope. Thus, in 
this embodiment of the present invention, the objective 
lens 20 is eliminated and replaced by the objective lens 
96 of the surgical microscope. The focal length of the 15 
lens 18 is adjusted to appropriately relay the Fourier 
plane 95 of the surgical microscope lens 96 to the aper- 
ture plane 30. In almost all other respects this embodi- 
ment is similar to the previously described embodiment. 
One possible exception is that if the F-number of the "® 
surgical microscope objective lens 96 is not sufficiently 
low to give the desired area of coverage on the cornea, 
then additional point sources 97 of light at multiple 
locations will be necessary outside the objective. These ^ 
additional sources may be created with an illuminated 
pinhole mask or optical fibers. 

FIGS. 11 through 14 illustrate an aspect of the system 
of the invention which enables both the epithelial cor- 
neal surface and the endothelial corneal surface to be 3Q 
detected and displayed in real time simultaneously. 
FIG. 11 shows an example of a reflected pattern 80 
which might occur at the detector 34 from the simple 
pattern shown in FIG. 6 comprising a pair of orthogo- 
nal linear arrays of light points. As indicated in FIG. 11 35 
each detected point 82 which is not on the optical axis 
will have a secondary reflection 84, of much lower 
intensity emanating from the back surface or endothe- 
lial surface of the cornea. The detected array might 
produce, for example, an intensity versus distance curve 4Q 
such as shown in FIG. 12. The long spikes 86 of light 
intensity represent the reflection of the discrete point 
light sources from the anterior, or front, surface of the 
cornea, with some degree of noise 88 occurring be- 
tween the spikes. A secondary spike or cluster 90 of 45 
light intensity which is discernibly higher than the noise 
88 occurs adjacent to each high intensity spike 86. This 
represents the lower-intensity reflection of the light 
points off the endothelial cornea surface. The plotting 
shown in FIG. 12 can easily be sampled or filtered to 50 
identify and separate the high intensity spikes 86 from 
the low intensity spikes 90. As can be appreciated by 
those skilled in this art, the programming in the com- 
puter can first determine the signal contribution from 
those spikes which achieve amplitudes above a prede- 55 
termined threshold and then subtracting the contribu- 
tion to the signal which correspond the high intensity 
spikes 86 to obtain a signal which contains only the low 
intensity spikes 90 and the noise 88. The process of 
identifying spike location for the high intensity spikes 86 60 
is now repeated for the low intensity spikes 90, but with 
a lower threshold. In some embodiments of the inven- 
tion, it may prove efficacious to electronically amplify 
the signal from which high intensity spikes 86 have 
previously been deducted in order to facilitate the 65 
threshold differentiation between low intensity spikes 
90 from the noise 88. It is important to note that this 
selection process is facilitated by the observation that 
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the specific amplitude of the spikes 86 and 90 are not as 
important as their actual location. 

FIG. 13 and 14 show separate plottings of intensity 
versus distance for the front surface, anterior reflection 
and the rear surface, posterior reflection. 

Once the primary and secondary reflections are 
known and located as in FIGS. 13 and 14. the shapes 
and elevation points of both the epithelial and endothe- 
lial surface can be calculated by the approximation 
method described previously, and two sets of data can 
thus be presented to the physician. Similarly, the cross 
sections and appropriate values can be represented in 
the lower two quadrants of the display as illustrated in 
FIG. 3. 

It should be understood that in the drawings and the 
description herein, as well as in the claims, references to 
"up", "down**, "Mower", "upper", "left" or "right" are 
intended only for convenience in referring to the em- 
bodiments as represented in the drawings, and not as 
limiting any possible orientations of the instrument or 
components. The drawing figures are not to scale. Fur- 
ther, the term "objective lens" as used herein and in the 
claims and drawing figures is intended to refer to either 
an objective lens specific to the instrument or an objec- 
tive or final focusing lens of a surgical microscope, if the 
instrument is used as part of a surgical microscope. 

The above described preferred embodiments are in- 
tended to illustrate the principles of the invention, but 
not to limit its scope. Other embodiments and variations 
to these preferred embodiments will be apparent to 
those skilled in the art and may be made without depart- 
ing from the spirit and scope of the invention as defined 
in the following claims. 

We claim: 

1. An ophthalmic diagnostic instrument for determin- 
ing the shape of the cornea of a patient's eye, compris- 
ing, 

an objective lens as an optical element of the instru- 
ment, on an optical axis of the instrument. 

means for projecting a pattern of discrete separated 
point light sources and forming a real image of the 
pattern of point light sources at a position located 
between the interior of the objective lens and the 
patient's eye, the real image including point light 
sources in positions which traverse substantially 
directly across the optical axis, 

means for expanding the region of coverage on the 
cornea by using the objective lens as a field lens for 
the pattern image, 

means for selecting and collecting a reflected image 
of the pattern as reflected paraxially off the cornea, 
and for detecting a reflected position of substan- 
tially each point light source, as reflected from the 
cornea, including means for relaying the Fourier 
plane of the objective lens to a relayed position in 
the instrument, with aperture means positioned at 
said relayed position for limiting the collected light 
to that which is reflected paraxially off the cornea, 
whereby the aperture means is a spaced distance 
from the objective lens, 

means for analyzing the returned, collected pattern 
image and for comparing it to the undistorted pat- 
tern as projected, including analyzing the relative 
location and spatial orientation of the reflected 
point light sources as compared to the pattern as 
projected, and 
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means for deriving mathematically a close approxi- 
mation of a corneal surface shape that would give 
rise to such a collected pattern image, 

whereby the real image point sources, extending sub- 
stantially across the optical axis, enable enhanced 5 
measurement of the central optical zone about the 
visual axis of the cornea. 

2. Apparatus according to claim 1. wherein the pat- 
tern of discrete separated point light sources includes 
more than one rectilinear sequence. 1° 

3. Apparatus according to claim 1. wherein the pat- 
tern of discrete separated point light sources includes an 
asymmetrical shape having a plurality of lines of point 
light sources. 

4. Apparatus according to claim 3, wherein the asym- 15 
metrical shape comprises generally a star with an odd 
number of points. 

5. Apparatus according to claim 1, wherein the means 
for expanding the region of coverage comprises means 
for placing the source image in the objective lens, thus 20 
using it as a field lens. 

6. Apparatus according to claim 1, in combination 
with a surgical microscope having a standard auxiliary 
camera mount, and the ophthalmic diagnostic instru- 75 
ment being connected to the surgical microscope via 
the auxiliary camera mount, with an objective lens of 
the surgical microscope serving as the objective lens of 
the ophthalmic diagnostic instrument. 

7. Apparatus according to claim 1, including an illu- ^ 
minating light source, a pattern plate or mask positioned 
for projecting light from the light source through the 
plate, beam splitter means in the path of projected light 
from the pattern plate for reflecting and folding the 
projected pattern into a path coaxial with the optical 35 
axis of the instrument, with the beam splitter means 
located in said spaced distance between the aperture 
means and the objective lens, and optical means be- 
tween the beam splitter means and the position of the 
patient's eye for focusing the projected pattern into the 40 
real image in front of the patient's eye. 

8. Apparatus according to claim 7, further including 
means for receiving a reflected pattern from the cornea 
back through the optical means and through the beam 
splitter means, detector means on the opposite side of 45 
the beam splitter means from said optical elements, 
further optical means for focusing the return reflected 
and distorted pattern onto the detector means, and said 
aperture being positioned in a Fourier plane of the re- 
turning reflected pattern to eliminate substantially all 50 
light from the detector means except light reflected off 
the cornea as parallel to the optical axis of the instru- 
ment, whereby the spatial orientation of the pattern 
detected on the detector may be compared to the origi- 
nally transmitted pattern for determination of the cor- 55 
neal shape through analysis of the positions of reflected 
points of the pattern. 

9. Apparatus according to claim 1, wherein the oph- 
thalmic diagnostic instrument includes means for fold- 
ing the pattern of discrete separated point light sources 60 
onto the optical axis of the instrument, toward the cor- 
nea, with the means for projecting the pattern including 

a source of the pattern off-axis from the optical axis and 
from the path of the returned, distorted pattern image. 

10. An ophthalmic diagnostic instrument for deter- 65 
mining the shape of the cornea, comprising, 

an objective lens as an optical element of the instru- 
ment, on an optical axis of the instrument, 



means for projecting a pattern of discrete separated 
point light sources and forming a real image of the 
pattern of point light sources at a position located 
between the interior of the objective lens and the 
eye. . . 

means for expanding the region of coverage on the 
cornea by using the objective lens as a field lens for 
the pattern image, 

means for selecting and collecting a reflected image 
of the pattern as reflected paraxially off the cornea, 
and for detecting a reflected position of substan- 
tially each point light source, as reflected from the 
cornea, including means for relaying the Fourier 
plane of the objective lens to a relayed position in 
the instrument, with aperture means positioned at 
said relayed position for limiting the collected light 
to that whiqh is reflected paraxially off the cornea, 
whereby the aperture means is a spaced distance 
from the objective lens. 

means for analyzing the returned, collected pattern 
image and for comparing it to the undistorted pat- 
tern as projected, including analyzing the relative 
location and spatial orientation of the reflected 
point light sources as compared to the pattern as 
projected, 

means for deriving mathematically a close approxi- 
mation of a corneal surface shape that would give 
rise to such a collected pattern image, and 

means for projecting a second light pattern compris- 
ing concentric circles toward the cornea simulta- 
neously with said pattern of discrete separated 
point light sources, and means for separately ana- 
lyzing distorted reflected light from the cornea 
relating to the concentric circles and for providing 
separate, qualitative information which can be 
compared with the corneal surface shape derived 
via the pattern of discrete separated point light 
sources. 

11. Apparatus according to claim 10, wherein said 
means for projecting a second light pattern includes a 
second light pattern illuminating source, a second pat- 
tern plate or mask, second beam splitter means posi- 
tioned along the axis of the instrument and in position to 
fold a projected pattern from the second pattern mask 
onto the optical axis of the instrument, and a first polar- 
izer in the path of the projected pattern of discrete 
separated point light sources and a second polarizer in 
the path of the second light pattern, with opposite po- 
larity established by the orientation of the two polariz- 
ers such that the two projected patterns are projected to 
real image locations with opposite polarity, and such 
that their reflections are more easily separable with a 
polarizing analyzer at the detector means. 

12. Apparatus according to claim 1, further including 
means for separately analyzing a secondary returned, 
reflected pattern image as reflected from the back or 
endothelial surface of the cornea. 

13. Apparatus according to claim 12, wherein said 
means for separately analyzing includes filtering means 
for electronically separating returned light points on the 
detector means occurring from the front surface of the 
cornea from those occurring from the back surface of 
the cornea, by separating different ranges of amplitude 
of the detected light. 

14. Apparatus according to claim 13, wherein said 
means for deriving mathematically includes computer 
means for determining the endothelial corneal surface 
shape from the locations on the detector means of the 
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detected light points reflected from the endothelial 
surface. 

15. Apparatus according to claim 1. wherein the 
means for projecting a pattern includes an illuminating 
light source and a plate with a laser-cut or photolitho- 5 
graphically produced pattern of discrete holes to form 
the discrete separated point light sources. 

16. An ophthalmic diagnostic instrument for deter- 
mining the shape of the cornea, comprising, 

an objective lens as an optical element of the instru- io 
ment, on an optical axis of the instrument, 

means for projecting a pattern of discrete separated 
point light sources and forming a real image of the 
pattern of point light sources at a position located 
between the interior of the objective lens and the 15 
eye. 

means for expanding the region of coverage on the 
cornea by using the objective lens as a field lens for 
the pattern image, 

means for selecting and collecting a reflected image 20 
of the pattern as reflected paraxially off the cornea, 
and for detecting a reflected position of substan- 
tially each point light source, as reflected from the 
cornea, including means for relaying the Fourier 
plane of the objective lens to a relayed position in 25 
the instrument, with aperture means positioned at 
said relayed position for limiting the collected light 
to that which is reflected paraxially off the cornea, 
whereby the aperture means is a spaced distance 
from the objective lens, 30 

means for analyzing the returned, collected pattern 
image and for comparing it to the undistorted pat- 
tern as projected, including analyzing the relative 
location and spatial orientation of the reflected 
point light sources as compared to the pattern as 35 
projected, 

means for deriving mathematically a close approxi- 
mation of a corneal surface shape that would give 
rise to such a collected pattern image, including 
means for utilizing the detected spatial orientation 40 
of the reflected point light sources to determine the 
corneal shape along a selected cutting plane on the 
eye in accordance with the general formula 

*^--(«^>[(«^) , + ']' 

whereby y is depth of a cornea reflection point 
from a datum plane, x is distance from the optical 50 
axis of the instrument, and (a, b) are the coordinates 
of the real image of an illumination point space. 

17. A method for determining the shape of the cornea 
of an eye, comprising, 

projecting a pattern of discrete separated point light 55 
sources and forming a real image of the pattern of 
point light sources at a position located in front of 
the eye, the real image including point light sources 
in positions which traverse substantially directly 
across the visual axis of the eye, 60 

selecting and collecting a reflected image of the pat- 
tern as reflected paraxially off the cornea, and de- 
tecting a reflected position of substantially each 
point light source, as reflected from the cornea, 

analyzing the returned, collected pattern image for 65 
comparing it to the undistorted pattern as pro- 
jected, including analyzing the relative location 
and spatial orientation of the reflected point light 
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sources as compared to the pattern as projected, 
and 

deriving mathematically a close approximation of a 
corneal surface shape that would give rise to such 
a collected pattern image. An ophthalmic diagnos- 
tic instrument for determining the shape of the 
cornea, comprising, 

an objective lens as an optical element of the instru- 
ment, on an optical axis of the instrument. 

means for projecting a pattern of discrete separated 
point light sources and forming a real image of the 
pattern of point light sources at a position located 
between the interior of the objective lens and the 
eye, 

means for expanding the region of coverage on the 
cornea by using the objective lens as a field lens for 
the pattern image, 

means for selecting and collecting a reflected image 
of the pattern as reflected paraxially off the cornea, 
and for detecting a reflected position of substan- 
tially each point light source, as reflected from the 
cornea, including means for relaying the Fourier 
plane of the objective lens to a relayed position in 
the instrument, with aperture means positioned at 
said relayed position for limiting the collected light 
to that which is reflected paraxially off the cornea, 
whereby the aperture means is a spaced distance 
from the objective lens. 

means for analyzing the returned, collected pattern 
image and for comparing it to the undistorted pat- 
tern as projected, including analyzing the relative 
location and spatial orientation of the reflected 
point light sources as compared to the pattern as 
projected, 

means for deriving mathematically a close approxi- 
mation of a corneal surface shape that would give 
rise to such a collected pattern image, and 

whereby the real image point sources, extending sub- 
stantially across the optical axis, enable enhanced 
measurement of the central optical zone about the 
visual axis of the cornea. 

18. The method of claim 17, wherein the pattern of 
discrete separated point light sources comprises a gen- 
erally cruciform shaped pattern with crossing rectilin- 
ear rays of point light sources as an intersection point 
lying on the optical axis of the instrument. 

19. The method of claim 17, wherein the pattern of 
discrete separated point light sources comprises a gen- 
erally asterisk shaped pattern with an intersection point 
at the optical axis of the instrument, and including 
means associated with the pattern for establishing a 
readily identifiable rotational orientation of the pattern. 

20. The method of claim 17, including using an off- 
axis illuminating light source and projecting light 
through a pattern mask and then reflecting the pro- 
jected pattern off a beam splitter to fold the projected 
pattern into a path coaxial with the optical axis of the 
instrument, focusing the projected pattern reflected off 
the beam splitter to form the real image in front of the 
patient's eye, receiving a reflected pattern from the 
cornea back through the beam splitter and focusing the 
reflected pattern onto the detector, and including pass- 
ing the returning reflected light pattern through an 
aperture en route to the detector to eliminate all light 
reflected off the cornea except that which is parallel to 
the optical axis of the instrument, whereby the spatial 
orientation of the pattern detected on the detector may 
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be compared to the originally transmitted pattern for 
determination of the corneal shape through analysis of 
the positions of reflected points of the pattern. 

21. A method for determining the shape of the cornea 
of an eye, comprising, 5 
projecting a pattern of discrete separated point light 
sources and forming a real image of the pattern of 
point light sources at a position located in front of 
the eye. 

selecting and collecting a reflected image of the pat- ^ 
tern as reflected paraxially off the cornea, and de- 
tecting a reflected position of substantially each 
point light source, as reflected from the cornea, 

analyzing the returned, collected pattern image for 
comparing it to the undistorted pattern as pro- 
jected, including analyzing the relative location 
and spatial orientation of the reflected point light 
sources as compared to the pattern as projected. 

deriving mathematically a close approximation of a 2Q 
corneal surface shape that would give rise to such 
a collected pattern image, and 

projecting a second light pattern comprising concen- 
tric circles toward the cornea simultaneously with 
the pattern of discrete separated point light 25 
sources, and separately analyzing reflected light 
from the cornea relating to concentric circles and 
providing separate, qualitative information which 
can be compared with the corneal surface shape 
derived via the pattern of discrete separated point 30 
light sources. 

22. The method of claim 17, further including sepa- 
rately analyzing a secondary returned, reflected pattern 
image as reflected from the back or endothelial surface 
of the cornea. 35 

23. The method of claim 22, wherein the step of sepa- 
rately analyzing includes electronically separating re- 
turned light point on the detector means occurring from 
the front surface of the cornea from those occurring 
from the back surface of the cornea, by separating dif- 40 
ferent ranges of amplitude of the detected light on the 
detector. 

24. The method of claim 23, including mathemati- 
cally deriving the endothelial corneal surface shape on 

a computer, from the locations on the detector of the 45 
detected light points reflected from the endothelial 
surface. 

25. A method for determining the shape of the cornea 
of an eye, comprising, 

projecting a pattern of discrete separated point light 50 
sources and forming a real image of the pattern of 



point light sources at a position located in front of 
the eye, 

selecting and collecting a reflected image of the pat- 
tern as reflected paraxially off the cornea, and de- 
tecting a reflected position of substantially each 
point light source, as reflected from the cornea, 

analyzing the returned, collected pattern image for 
comparing it to the undistorted pattern as pro- 
jected, including analyzing the relative location 
and spatial orientation of the reflected point light 
sources as compared to the pattern as projected, 
and 

deriving mathematically a close approximation of a 
corneal surface shape that would give rise to such 
a collected pattern image, 
including utilizing the detected spatial orientation of 
the reflected point light sources to determine the 
corneal shape along a selected cutting plane on the 
eye in accordance with the general formula 

where y is depth of a cornea reflection point from 
a datum plane, x is distance from the optical axis of 
the instrument, and (a, b) is the coordinate location 
of the real image of an illumination point in space. 

26. The method of claim 17, further including elec- 
tronically producing cross sectional images of the eye at 
selected cutting planes using information derived from 
the mathematical derivation of the corneal surface 
shape. 

27. The method of claim 17, wherein the pattern of 
discrete separated point light sources is projected 
through the objective lens of a surgical microscope, and 
including collecting the returned pattern image through 
the objective lens of the surgical microscope. 

28. The method of claim 21, further including pro- 
ducing and displaying an image showing the distortion 
of the projected concentric circles. 

29. The method of claim 17, wherein a front objective 
lens comprises an element closest to the eye of the pa- 
tient, and including spacing the front objective lens at 
least about 110 mm away from the eye. 

30. The method according to claim 17, including 

forming the real image of the pattern of point light 

sources substantially in the objective lens, so as to use 

the objective lens as a field lens, enabling a larger field 

of view of the cornea. 
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ABSTRACT 



The present invention describes a technique and appara- 
tus for finding spots in an image with substantial noise 
making it difficult to identify without specialized noise 
suppression algorithms. In the context of determining 
corneal shape, as an example of the technique, the re- 
flections of point light sources in or on the cornea have 
long played a diagnostic role. The image analysis tech- 
nique described applies the tools of mathematical mor- 
phology and prior information about the shape of illum- 
ination patterns to remove noise and isolate the points of 
.interest for further mathematical analysis. The output 
from the technique is a set of pairs matching the de- 
tected points in the image with the known location of 
the illumination. 

17 Claims, 7 Drawing Sheets 
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The process of separating the spots of reflected signal 
APPARATUS AND METHOD OF IDENTIFYING from the background light and locating the peak or 
SIGNALS IN BIOLOGICAL TISSUES center of the signal presents a challenging image pro- 
cessing problem which can be crucial to the success of 
REFERENCE TO RELATED APPLICATION 5 the entire surface profiling procedure. The apparatus 
This application is a continuation-in-part of copend- and methods for obtaining the accurate location of such 
ing U.S. application Ser. No. 456,109, filed Dec. 22, S1 8 nals are the sub J ect matter of this invention. These 
1989, now U.S. Pat. No. 5,054,907. techniques are by no means restricted to ophthalmic 
" " applications, or even biological systems. These types of 
BACKGROUND OF THE INVENTION 10 techniques have been in use in military applications 
The invention relates generally to optics, and more v y he [ e t f r S ets are at times difficult to differentiate from 
particularly to an imaging system wherein noise filter- the background field and where filtering techniques 
ing techniques are employed to increase the signal to such as thresh olding and fast Fourier transform filtering 
noise ratio used in analyzing the image. P rove ""successful. medicine and, m particular, m 
In many biological systems, a method of preference 15 surgery ' co fP uters h f ave only recently begun to be 
for identifying the location, shape, or other characteris- lnc ° r PO r ated as part of surgical devices. Since the tech- 
tics of a tissue is to shine light or some other form of nl £ u f S t that a 5 e , th , e suh £ ct of lnv£ f ,on 
radiation onto the tissue and observe either how the sub * an * lal calculation, these techniques have only re- 
light is reflected, absorbed, refracted, or scattered by , n ^a^Th^L uY envir ° nme f , 
the tissue. Inherent to observing the radiative effect on 20 .f ,^ furtber be l™.' the a Pf aratus and 

said tissue is the ability to detect a signal emanating °t ^ *T i T ^ ^ 

a.™ „ ^ -r a' + u- ^ filtering techniques of mathematical morphology many 

trom the tissue as a consequence of disturbing the tissue r t. - u ■ * j i / T *r 

• „ , ■ . t /• i • -\m , . . 6 4 , of which onginated with the pioneering work of J. Von 

m some designated fashion. Many biological structures, x T ~ r j i - * * j j • * i 

. Z * v *l r r Neumann for developing automated devices to analyze 

however, are very sensitive to light or other forms of • t • • i ■ . 

,. v. J\, r . , ;T 25 images by companng a given pixel of the image with its 

radiation thus limiting the threshold of signal that can imm W e neig W>rs. Ma ny of these nonlinear filtering 

be used to generate nondamaging responses. Even more techni s are described b Serra (J Aml is J d 

troublesome is that many biological tissues have weak Mathematieal Morphology, Jean Serra, Academic Press 

reflectivity or weak scattering properties so that the 19g2) and borrow frQm the fidds of al braic topo logy, 

pertinent signal resulting from the disturbance is diffi- 30 harni onic analysis, stochastic processe^ntegralgeome- 

cult to detect. This presents the observer with the com- t and Qthers We , y these nonlinear fllteri tech . 

pound problem of having to detect weak signals and ni of mathematical morphology to the problem of 

lo ^ na ! 1 /° t nolse / 1 f tlos - u j* • , . Elating and identifying the centers of refiected point 

This is illustrated herein by discussion of the anterior light sources In the case of the invention deS cribed by 

Tu u t, . eyeS C ° rnea i kn ° Wn aS « he e P Ithelmm > 35 Sklar et al. (U.S. patent application Ser. No. 456,109, 

although this discussion applies to any reflecting surface now tj. s . Pat . No. 5,054,907) the purpose for identify- 

withm the eye. When trying to define the surface topog- ing said re fl e cted light sources, herein after referred as 

raphy of the epithelium, it is often desirable to shine "p oin ts", is to present the data to a surface profiling 

pomt sources of light onto the epithelium from a pre- a i gori thm to determine the topography of the corneal 

cisely established location and to measure very accu- 40 surface . l n other uses of the present invention not only 

rately on a detector the location of the reflected image different structures are to be described, but different 

from said point source. However, point sources are a types Q f image sources and scattering as well as absorp- 

distnbution of energy about some given location and tionSj reflections, or refractions can be considered, 

said distribution is not generally constant or uniform The problem of identifying the location of the centers 

from source to source. Since these energy distributions 45 Q f the points is addressed in three stages. First, the 

are reflected from unspecified surfaces at unpredeter- po j nts must be isolated. Then we must look for and 

mined angles of incidence and with varying reflectivity recognize symmetry in the pattern of points. And fi- 

and scattering, the identification of the reflected images na n y , we have rejection of noise pulses which have 

of the point sources can be a difficult problem, espe- eluded previous elimination criteria, 
cially when low illumination levels are desired. 50 

The accuracy of the technique used by Sklar et al. to SUMMARY OF THE INVENTION 
measure the surface profile of the cornea (U.S. patent in accordance with the present invention, an appara- 
application Ser. No. 456,109, now U.S. Pat. No. tus and method for identifying the location of the cen- 
5,054,907) hinges on the ability of the apparatus to mea- ters of signals with poor signal to noise ratios and poor 
sure the location of the rays from point sources in front 55 background differentiation are described. The method 
of the eye as they are reflected from the cornea. As calls for first isolating the signal "points", recognizing 
described by Sklar et al., the reflected rays are detected the symmetries in the point patterns, and suppressing 
by an imaging device such as a CCD camera and digi- noise signals which have eluded previous elimination, 
tized using a frame grabber card in a computer. Other Since the recorded image background is frequently 
perimeter devices used to measure eye features also rely 60 nonuniform due to scattering from the instrument, dust 
upon accurately measuring the response of the eye to particles in the air, or even from the target, the first step 
shining light onto the eye. In general, the observed in isolating the point signals is subtraction of the back- 
reflections of these spots of light are embedded within a ground. This is accomplished by subtracting a morpho- 
nonuniform background which may be very close to the ~* logical gray-scale opening of the image using a circular 
noise floor level for contrast selectivity or threshold 65 kernel, which corresponds to an erosion followed by a 
filtering because of poor light source edge definition dilation, from the original picture. The terms "erosion", 
and because of low reflectivity from the corneal epithe- "dilation", "opening", "closing", "morphological oper- 
lium. ations", "circular kernel", "hit or miss topology" and so 
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forth are intended in their algebraic sense as defined by onto the 12 millimeters which defines the diameter of 

J. Serra as cited above. For example, given two finite the cornea taking into account the magnification of the 

sequences Gi, . . . , G/, . . . , G m of open sets and Ki, . apparatus as seen by the detector. In practice, the kernel 

. . Kfc, . . . , Kji of compact sets in a real space of N radius should not exceed approximately half the aver- 

dimensions, the class of all closed sets which hit every 5 age separation between intensity peaks. We say average 

Gi and miss every K/, defines an open neighborhood in because in general a uniform distribution of the peaks is 

the space F of all closed sets in a real N dimensional neither required or optimally desirable, 

space. For reasons of computational speed, it is preferable to 

In the present invention, the size of the kernel is make the kernel radius as small as possible. One of the 
chosen so as to remove all features smaller than the 10 limitations to reducing the kernel radius is that if the 
approximate size of the dots from the image. The size of peaks are too jagged, the difference operation described 
the dots is a known quantity since it is related to the above will reveal a variety of local minima rather than 
illumination source. In the embodiment of the invention a single point at the peak of the spot whenever the 
described below this corresponds to using a circular smaller peaks are separated by more than the kernel 
kernel of radius 3 pixels since the point sources are 15 radius. This proliferation of intensity peaks leads to 
being magnified by the internal optics of the system to confusion and degradation of the result. Fortunately, 
reflections which appear of that size to the detector we can test for adequate choices for the kernel radius 
observing the reflections. This first level of filtering of during experimental phases of setting up the apparatus 
the original image, hereinafter referred to as Io, shown for a given application in order to avoid the problem of 
in FIG. 1, leaves the resultant image, hereinafter re- 20 proliferating peaks. Alternatively, the apparatus is set 
ferred to as Ii, shown in FIG. 2, with a level back- up to allow for modification of the kernel radius when- 
ground and tends to equalize the maximum light intensi- ever the resulting image definition is insufficient be- 
ties observed from the dots by reducing the more in- cause of the detection of excessive number of intensity 
tense spots proportionately more than the fainter spots. peaks. For corneal topography applications we use a 

The next step in the process of isolating the points is 25 kernel of radius 5 pixels, 

to locate the peak intensities of the spots. This operation The resulting image I2 described above as a result of 

involves a combination of dilation with thresholding. A the morphological transformations is to produce a list of 

copy of the result of the first level of filtering described points in l2made up primarily of the light rays originat- 

above, I], is stored in computer memory while a second ing in Ioand reflected from the surface of the cornea. It 

image, I2, is generated by applying a threshold transfor- 30 is anticipated that I2 will also contain some spurious 

mation to 1 1. As a consequence of this threshold filter- noise spikes. 

ing transformation, all pixels that composed Ij which The next step in the process is to use knowledge of 

were below the threshold value are set to the threshold light source distribution symmetries to further suppress 

value while those above the threshold are left un- noise in I2. The operation now is to find the center of an 

changed. 35 approximately symmetrical pattern so that it can be 

The threshold value is determined experimentally used to segregate a number of lines that pass through 
and is determined according to the types of light this center into separate groups. We say approximately 
sources, optical elements, and targets to be identified. symmetrical patterns for two principal reasons. First, 
The threshold value is chosen to be low enough to be the optics of the system and the reflecting surface itself 
below the expected intensity level of the spots originat- 40 will introduce some level of pattern distortion into I2. 
ing from the source light pattern, but above most of the Second, there is no requirement that we restrict discus- 
observed random noise emanating from the floor. A sion to only one pattern or even one type of pattern, 
dilation is then applied to the image I2 once again using Thus, to simplify the analysis a sorting into distinct 
a circular kernel. The effect of this second filtering pattern groups is desirable. An example of how these 
operation is to leave the maximum intensity of each spot 45 points might be distributed is shown in FIG. Sa. 
unchanged and at the same location while spreading In one embodiment of the present invention, the oph- 
said maximum value of the spot around the peak loca- thalmic application described herein, a star shaped dis- 
tion outwards to the full extent of the radius of the tribution of light sources being projected through the 
circular kernel. These peak values are determined by ' objective lens of a surgical workstation (Sklar et al, U.S. 
subtracting Ii from I2 and looking for areas in the result- 50 patent application Ser. Nos. 475,657 and 456,109, now 
ing image whose intensities are identically zero. These U.S. Pat. No. 5,054,907) was selected as the desired 
areas of zero intensity correspond to the regions in mask in order to provide high definition of corneal 
image I2 that were above threshold and represented curvatures in the neighborhood of the visual axis of the 
local intensity maxima within the region determined by eye as well as a measure of angular asymmetry, 
the radius kernel since they were unchanged between 55 Even though the light source or input template used 
the original image Ioand the dilated image. The effect of to produce the illumination pattern is symmetrical in 
these operations is shown in FIGS. 3 and 4. this embodiment of the invention, the output pattern as 

The choice of kernel radius for the dilation kernel observed by the detector is typically not symmetrical, 

must be based on the separation distance between adja- Variations in illumination intensity which give rise to 

cent spots in the light source which gave rise to Io and 60 asymmetries can be caused by difference in fiber cou- 

the anticipated smoothness of each of the intensity pling the light source to the fiberoptic waveguides used 

peaks observed. In practice, an illumination mask or a to establish a pattern, or differential losses within the 

pattern of light sources is used with a given separation waveguide, or pinhole size variations is fabricating a 

distance between the originating light spots. In ophthal- mask for the illumination source. Such variations can 

mic applications where you are investigating the shape 65 lead to loss of detectable spots in the image I2. Other 

of the cornea, the distance between light spots is se- mechanisms for spot loss can be ascribed to variations in 

lected to correspond to the desired resolution as pro- reflectivity from the target being examined and even 

jected through the optical instrument of the apparatus from image processing. The noise spots that have either 
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been created through the image processing or which shown in FIG. 6. We thus divide the spots observed in 
have survived the various filtering operations described I2 into eight separate bins using angles midway between 
above will not generally be distributed symmetrically. each known line as the bin boundaries. So far, we have 
Even more important, the distribution of noise spots are used the knowledge that the illumination light sources 
unlikely to correspond to the preselected source sym- 5 have been arranged along discrete and distinct patterns 
metries even when distorted. This information is used in and we have used the existence of such patterns to 
further suppressing noise spots. locate the center of the pattern and to sort the observed 

The process for finding the center of the pattern spots into distinct bins or subpatterns. We now proceed 
begins by going through the pattern and assuming that to use knowledge of the original shape of the distinct 
each peak in the pattern is the center of the pattern. The 10 illumination pattern to further suppress noise, 
assumption that each peak is the center of the pattern is The final noise rejection operation in the system is 
then tested by examining the reflection of each other based upon a priori knowledge of the configuration of 
spot in the pattern about the supposed center and the illumination system. In the example noted above 
searching for existence of another spot close to the concerning topographical mappings of the eye's cornea, 
reflection. Namely, the vector from the assumed center 15 the subpatterns described in FIGS. 9 and 11 correspond 
to the test point is reversed. If there exists another point to straight lines intersecting at some center point. The 
near the reflection, the given score scale value for that angles of each line with respect to one another are pre- 
particular assumed center is increased, After sequen- determined during fabrication of the illumination masks, 
tially going through this process of testing each peak for Using this information, the spots in image I2 were segre- 
closest neighbors with every surviving spot in I2, the 20 gated into individual groups or bins. Each grouping of 
peak with the highest score corresponds to the center. points in a bin can then be approximately rotated to the 
This process is illustrated in FIG. 5b. horizontal axis by making use of the knowledge of the 

This process of locating the center is calculationally angular relation between each line in the illumination 
time consuming. It is, however, very robust. It success- mask. 

fully locates the center of a given illumination pattern as 25 Any point within the group lying more than a prese- 
reflected by the eye even when many spots are missing lected number of pixels off the horizontal axis following 
or with a multiplicity of extraneous noise spots. When the rotation is identified as a noise pixel and deleted 
dealing with the variety of ailments which are often from the image. This operation is shown in FIG. 6. In 
encountered by ophthalmologists, keratoconus being practice, we have found that a discrimination value of 5 
one notable case, surface aberrations can frequently 30 pixels is adequate for the task of noise suppression. After 
cause illumination spots to go astray. A robust system of elimination of such noise pixels, the subpatterns can be 
detection which will allow an instrument to perform counterrotated to its corresponding alignment to form 
irrespective of such variation was one of the motiva- image I3. In general, there is no requirement that the 
tions in the design of the present invention. With the subpatterns correspond to straight lines or that the su- 
computational capacity available today with high per- 35 perposition transformation be isometric. It is important 
formance dedicated computer boards using 80386 and however that the patterns be selected so as to permit 
80486 based microprocessor with specialized mathemat- superposition transformations which are invertible. 
ical coprocessor chips and vectored accelerator boards, Once the points that are well away from the true line 
the calculational time of the processes described herein have been eliminated, the resulting pattern I3 is com- 
which were formerly practicable on mainframe com- 40 pared with a stored calibration pattern where the scales 
puters become accessible to commercial uses and surgi- of the two patterns are matched by comparing the aver- 
cal environments. age spacing between points in each line subpattern. 

The template pattern 46 described in FIG. 9 and the Points which do not have a mate in the calibration 
pattern 82 in FIG. 11 of Sklar et al. (U.S. application pattern are further ignored before proceeding to the 
Ser. No. 456,109) correspond to a sequence of light 45 profiling procedure. Examples of such profiling proce- 
sources arranged along a number of intersecting dures are discussed in Sklar et al U.S. Pat. No. 
straight lines. If we apply the above general discussion 5,054,907. 

to these templates, the effort in finding the center of the _____ T __ T __ _ ____, „___ T _ 0 

pattern reduces to locating the intersection point for the ^ DESCRIPTION OF THE DRAWINGS 
straight line patterns as observed in image I2. Once the 50 FIG. la shows the illumination pattern Ioof the light 
center of the pattern is established, the approximately reflected from the cornea or from the tear layer overly- 
straight lines of spots radiating outward from this center ing the cornea as digitized by the frame grabber. A trace 
in image I2 are easily identified by collecting spots to- across the pattern at 1 shown in FIG. lb illustrates the 
gether that have approximately the same angle relative nonuniform peak heights and background, 
to the axis of the line based upon the arctangent value at 55 FIG. 2a shows the illumination pattern Ii after the 
the given spot. background has been subtracted. FIG. 2b shows a trace 

The grouping of spots into independent straight lines across th* pattern illustrating the equalization of the 
can be based upon the a priori knowledge of the number background level. 

of such straight lines which composed the initial illumi- FIG. 3a shows a trace 3 of image Ii, a reasonable 
nation pattern as described in FIG. 9 and 11 noted 60 threshold value 2, and a trace 3a across the dilated 
above. Thus, the correct number of angular bins can be threshold image. The difference between trace 3 and 
chosen for the collection system as' well as reasonable trace 3a is shown in FIG. 3b with trace 4 representing 
values for the boundaries of the bins. In other words, we the threshold level and the points 5 representing the 
are using knowledge of the shape of the illumination points at which the difference is identically zero indicat- 
pattern in order to improve noise suppression in a sys- 65 ing the local maxima of trace 3. 

tern with poor signal to noise ratios. FIG. 4 shows the resulting image after the difference 

In our illustrative example, we use eight straight lines operation with the threshold set so that only the null 
of illumination source points passing through the center points are visible. 
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FIG. 5a shows the distribution of points found in 
FIG. 4. The points labeled 2 represent noise signals. 
FIG. 5b illustrates the process of finding the point about 
which the pattern is symmetrical. A trial center point is 
denoted by 56. The vector 5 denotes the vector from 5 
the trial center point 5b to a given test point 5c. The 
vector 6a depicts the reflection of vector 6 about the 
trial center 5b and the point 5d is a point near the re- 
flected vector 6a, thus incrementing the score of point 
5b. 10 

FIG. 6a shows the pattern of dots with angular bins 
as indicated in trace lc. The horizontal axis is 7. To 
remove noise spots from one of the lines, say line la, the 
pattern lc is rotated through angle 8 to give the result 
shown in FIG. 6b. The lines 9 show the discrimination 15 
threshold, in our case 5 pixels, and the points 5c repre- 
sent points that are rejected since they lie sufficiently far 
off line, 

FIG. 7 shows an embodiment of an apparatus which 
illuminates the target and provides the optical means for 20 
detecting the signal and sending it the microprocessor 
means that initiates the noise suppression techniques 
discussed herein. 

FIG. 8 is a schematic drawing of the information flow 
for the filtering techniques that are a part of the present 25 
invention. The frame grabber board 101 inputs the 
image Io to the image opening operation 102. We then 
develop at 103 the image of Io with the background 
removed. The thresholding operation 104 replaces all 
intensity values less than the threshold by the threshold 30 
value and is followed by the dilation operation 105. The 
pattern 106 selects the points that are identically zero 
representing the local maxima in Io. The process of 
finding the center of the pattern 107 uses symmetry 
scoring. Next follows the process 108 of segregating the 35 
spots into lines based on the a priori information about 
the angles of the lines in the illumination source. The 
process 109 of rejecting the spots which are off the line 
after rotation then follows. The final step 10 denotes the 
process of matching the spots in the image lines to spots 40 
on the target lines. Non-matching spots are deleted. 

DESCRIPTION OF PREFERRED 
EMBODIMENTS 

In the drawings, FIG. 7 shows in schematic represen- 45 
tation an example of how an instrument would make use 
of the noise suppression features and methods of the 
present invention. It describes a system of optical ele- 
ments in accordance with the invention for use in carry- 
ing out ophthalmic diagnosis and analysis similar to a 50 
device taught by Sklar et al. (U.S. patent application 
Ser. No. 456,109). 

The system, generally identified by the reference 
number 110, includes an illuminator or light source 12, 
a pattern plate or disk 14 having a pattern of holes cut 55 
in the plate for producing a desired pattern of discrete 
light sources, a non-distorting beam splitter 16, a lens 18 
which projects an image of target 14 onto an image 
plane at 22. This image plane 22 is close to or coinciden- 
tal with the system of objective lens 20. The purpose of 60 
placing the image at this location 22 is to have the ob- 
jective lens 20 serve as a field lens, that is bending the 
rays of light 22a that form the image towards the pa- 
tient's cornea 21. 

As indicated in FIG. 7, the focused image 22 of the 65 
pattern is a real image, formed at some plane at or near 
the lens 20 and between the lens and the patient. The 
real image preferably is in the lens 20, but it can be very 
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closely in front of the lens (i.e. a few millimeters in 
front) or behind the lens. The importance is that the 
focused image 22 lie along the optical path passing 
through the lens 20 which can give optimal definition 
along the visual axis of the cornea. In this real image, 
each point source of light 22a projects a cone of light 
towards the patient. Thus, each point source 22a in the 
real image makes an infinite number of specular reflec- 
tions off the front surface of the cornea 24 of a patient's 
eye 26. In addition to the focused image 22, it may 
prove desirable to extend the range of coverage of a 
curved reflecting surface 21 by projecting a longer 
array of image points 22a than can be projected through 
the final lens 20. 

Towards this end, a light source 41 illuminates a 
fiberoptic bundle 42 whose polished terminations 43 
serve as additional image points to the image points 22a 
located within or directly in front of the final lens 20. 
The purpose of light points 43 is to increase the length 
of the image line 22 and thus provide more extensive 
peripheral definition of cornea 21 away from the visual 
axis of the eye. The position of light points 43 is adjusted 
so as to lie along the plane of pattern 22. 

Alternatively, even for conventional perimeter de- 
vices that only utilize light sources as described by light 
points 43, the present invention provides a means for 
improving the detection of said light points 43 when 
poor signal to noise ratios are perceived by the detector 
34. This enables weak light sources 12 and 62 to be 
utilized without sacrificing resolution, and thus mini- 
mize discomfort to the patient or light toxicity to the 
patients eye. 

The F-number of the final lens 20 determines the 
maximum area of the cornea that can be measured. The 
objective lens serves as a field lens, and the patient's 
cornea must be at the focal length of the lens 20. A 
preferred method for assuring that the eye is located at 
the focal length of the lens 20 is discussed in Fountain et 
al. (U.S. patent application Ser. No. 655,919). This as- 
sures that the light reflected off the eye parallel to the 
optical axis of the instrument is then brought to a point 
behind the lens 20 at the focal distance of the lens 20. 
This enables the return light to be apertured down, to 
select only those rays which were paraxial off the eye. 
This enables the system to localize a detected point to a 
point on the cornea from which that ray is reflected. If 
the objective lens 20 were not situated to serve as a field 
lens, outermost points of light in the pattern would not 
reflect off the cornea. As a field lens, the lens 20 effi- 
ciently bends the outer points of light toward the eye. 

It is preferred that the focal length of the lens be great 
enough to provide an unobstructed, comfortable dis- 
tance from the instrument to the patient and adequate 
working room for the surgeon, for surgical applications. 

FIG. 7 schematically indicate-s that the detector or 
camera plane 34 is connected to a microprocessor 50 
which may contain a frame grabber board 51 which 
digitizes the images detected by detector 34 together 
with a vector processor 52 for speeding the calcula- 
tions. The microprocessor may be connected to a dis- 
play device, such as a CRT monitor 53 as indicated. 
Data gathered from the system as described is received 
by the microprocessor 50 and analyzed. As described 
below, each detected point is correlated with the loca- 
tion of the particular point in the source pattern from 
which it emanated. This information is then provided to 
an algorithm as described by Sklar et al. (U.S. patent 
application Ser. No. 456,109, now U.S. Pat. No. 
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5,054,907) for display as part of a user interface, for for symmetry involves reflecting the vector 6 about the 

ophthalmic surgical procedures such as described by point 5b and looking to see whether there is another 

Sklar et al. (U.S. patent application Ser. No. 307,315 point in the pattern near this reflection of the vector. In 

now U.S. Pat. No. 5,098,241 and Ser. No. 475,657) or the case depicted, there is indeed a point nearby, namely 

for diagnostic verification of the eye's topography. 5 point 5d, so the symmetry score for point 5b would be 

FIG. 8 shows a chart of the flow of information and incremented. After using all the remaining points in the 
the multiple steps taken in filtering the noise from an pattern as test points for the assumed center 5b, we 
image such as those contemplated in the present inven- would arrive at a total symmetry score for the assumed 
tion. The chart, discussed above, traces the various center 5b. After checking each point in the pattern as an 
mathematical morphological transformations which are 10 assumed center, we arrive at a symmetry score for each 
applied to the images generated in the particular case of point. The point with the highest symmetry score is the 
the configuration shown in FIG. 7 which constitutes an actual center of the pattern. While the process of find- 
application in the field of ophthalmology. ing the center of the pattern is a very simple operation 

As explained above, FIG. 8 describes the flow of for the human eye, even when points are missing and 

information from the frame grabber 101, noted as the 15 extraneous points find their way into the field of view, 

frame grabber means 51 in FIG. 7, to the profiling sys- this is a nontrivial process for a computer. The process 

tern. The image from the frame grabber 101 is shown in. described as part of the present invention is robust even 

FIG. la. It is formed by the light from the spots 22 in in the presence of noise. 

FIG. 7 reflecting from the cornea or the cornea's tear Once the center of the pattern has been identified, the 

layer, passing through the aperture, and hitting the 20 angles of the pattern lines relative to the axis, known 

detector. A trace across this pattern along the line from the construction of the target, noted as 14 in FIG. 

shown as 1 in FIG. la shows the intensity along the line 7, can be used to segregate the dots into distinct lines, 

as a function of position. The result is shown in FIG. 16. This process is depicted in 108 of FIG. 8. FIG. 6a shows 

The pattern in FIG. lb is referred to as the image Io. the pattern with segregated lines lc drawn at the bisec- 

A copy of Io is passed through the morphological 25 tors of the pattern lines. Any points that lie between 

opening operation using a circular kernel shown as 112 adjacent segregation lines are considered to belong to 

in FIG. 8. The difference between the original image the same line. The set of points 7 form the line that is on 

land this opening is represented by 103 in FIG. 8. We the horizontal axis of the system, while a line la corre- 

refer to this image as Ij. The effect of the opening and spends to a line off the horizontal axis at an angle ft, 

the difference is to equalize the background level 30 shown as 8. Note that ft is known from the fabrication 

throughout the image, thus effectively removing the of the target illuminator. 

background from the resulting image shown in FIG. 2a Each group of dots grouped as a line by process 108 

with a trace across the pattern illustrating the equaliza- of FIG. 8 is examined for extraneous dots. This is done 

tion of the background level when compared with FIG. in 109 of FIG. 8 by rotating each group through its 

lb- 35 corresponding angle ft and deleting points that are 

While a copy of Ij is maintained, the information is more the discrimination value away from the horizontal 

passed first through a thresholding operation in which axis. The discrimination levels 5 are used to keep points 

all pixel values in the image less than the threshold within the lines and delete points 6 which fall outside 

value are set to the threshold value, represented sche- the criteria. This process is performed for each line in 

matically as 104 in FIG. 8, the output of which passes 40 turn deleting points outside a given distance from the 

through a dilation operation 105 in FIG. 8. The effect of axis. 

these operations is shown in FIG. 3a. In FIG. 3a, the The final step in the process involves matching the 
line 3 shows a trace across Ij, and 2 represents a reason- detected points in the image with the known location of 
able threshold level. Trace 3a shows what the trace points in the illumination target 14 of FIG. 7. This pro- 
looks like after both thresholding and dilating. The 45 cess is represented by 110 in FIG. 8. Since the center 
difference between \\ and the thresholded and dilated point is known from process 107 of FIG. 8, it is the first 
image is noted in 106 of FIG. 8. This difference is re- point to be matched. Since the known angle of each of 
ferred to as image I2 and is shown in FIG. 4. The places the target lines was used to get rid of noise pulses during 
in hthat are identically zero represent the local maxima" the rotation of image lines, it can also be used to tell 
of Io, which are the points sought by this technique. A 50 which group of dots goes with which line of the illumi- 
trace across this pattern is shown in FIG. 3b with the nation target 14. While the magnification between the 
points labeled 5 being the places that are identically image and the target is unknown since it depends on the 
zero in the pattern corresponding to the local maxima of radius of curvature of the cornea, the average spacing 
3 in FIG. 3a. of the dots on the illumination target line and the dots 
The collection of points in I2 that are zero must then 55 on the image line can be used to scale each line so that 
be organized for subsequent analysis. These points are the two patterns are normalized. Thus the distances of 
shown in FIG. 5a. Because of noise in the original input the points :on both the target and the image from the 
image, there are still a few noise spots in the pattern center are divided by their respective averages. The 
illustrated by 5a in this figure. These points are elimi- normalized points from both patterns are then examined 
nated by the following procedure. The first step in the 60 to make sure that there is a corresponding point within 
process is to find the center of the pattern represented a sufficiently small range chosen to be small compared 
by 107 in FIG. 8. The process assumes that each of the to the normalized dot spacing. Any point in either the 
spots I2 is the center of the pattern in turn and then image or the target pattern that is found to be without a 
examines the rest of the spots in the pattern to measure corresponding mate is discarded before the matched 
if they are symmetrical about the assumed center as 65 data set is sent to the profiling algorithm described by 
illustrated in FIG. 5b. Point 5* denotes the assumed Sklar et al. (U.S. Pat. No. 5,054,907). 
center. Point 5c is a test point. Vector 6 is the vector The above described preferred embodiments are in- 
joining point 5c to point 5c. The operation of checking tended to illustrate the principles of the invention, but 
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not to limit its scope. Other embodiments and variations 7. Apparatus according to claim 1, further including 

to these preferred embodiments will be apparent to means for projecting a variety of predetermined light 

those skilled in the art and may be made without depart- patterns whose reflected image can be individually 

ing from the spirit and scope of the invention as defined mapped isometrically onto a straight line simulta- 
in the following claims. 5 neously with said pattern of discrete separated point 

We claim: light sources, and means for separately analyzing dis- 

1. An ophthalmic diagnostic instrument for determin- torted reflected light from the cornea relating to the 
ing the shape of the cornea, comprising, invertible shapes and for providing separate, qualitative 

an objective lens as an optical element of the instru- information which can be compared with the corneal 
ment, on or symmetrically about an optical axis of 10 surface shape derived via the pattern of discrete sepa- 
the instrument, rated point light sources. 

means for projecting a pattern of discrete separated 8 - Apparatus according to claim 1, further including 
point light sources and forming a real image of the means for separately analyzing a secondary returned, 
pattern of point light sources at a position located reflected pattern image as reflected from the back or 
between the interior of the objective lens and the 15 endothelial surface of the cornea. 
ey C) 9. Apparatus according to claim 1, further including 

means' for expanding the region of coverage on the means for separately analyzing a secondary returned, 
cornea by using the objective lens as a field lens for reflected pattern image as reflected from the anterior 
the pattern image, surface of the e y e ' s lens - 

means for selecting and collecting a reflected image 20 10 Apparatus according to claim 1, further including 
of the pattern as reflected paraxially off the cornea, m * ans for se P a rately analyzing a secondary returned, 
and for detecting a reflected position of substan- reflected pattern image as reflected from the posterior 
tially each point light source, as reflected from the sur / ac e of the eye s lens. 

cornea, including means for relaying the Fourier „ 11 A PP aratus according to claim 1, further including 
plane of the objective lens to a relayed position in 25 m f ns *> r se P ara ! el V analyzing a secondary returned, 
the instrument, with aperture means positioned at re ^ Ct f d P attern iraa ^ as reflected from the retina 
said relayed position for limiting the collected light 12 A /P aratus according to claim 1, further including 
to that which is reflected paraxially off the cornea, m * anS f f ? T s fP aratel y analyzing a secondary returned, 
whereby the aperture means is a Spaced distance , Q "J^J ^e™^ * " 

from the objective lens, A . A t . 0 , ., 

f *i * , „ 4 , 13. Apparatus according to claim 8, wherem said 

means for analyzing the returned, collected pattern me ans for separately analyLg includes filtering means 

image and for filtenng the noise from the pattern for electron £ ally separating returned light points on the 

image using mathematical morphological transfor- detector means occurring from the fro * t s ; rface of the 

mations, 35 cornea f rom t jj Qse occurr i n g f rom the back surface of 

means for comparing the filtered, collected pattern the C0Tnea> b separating different ranges of amp i itude 

image to the undistorted pattern as projected, in- of the dete cted ij ght 

eluding analyzing the relative location and spatial u Apparatus according to claim 1, wherein the 

orientation of the reflected point light sources as means for pro j ecting a pattern i ne ludes an illuminating 

compared to the pattern as projected, and 40 Hght source and a plate with a laser „ cut or pno tolitho- 

means for deriving mathematically a close approxi- graphically produced pattern of discrete holes to form 

mation of a corneal surface shape that would give the d i scre t e separated point light sources, 

rise to such a collected pattern image. 15 A method for determining the shape of the cornea 

2. Apparatus according to claim 1, wherein the pat- 0 f an eye , comprising 

tern of discrete separated point light sources includes 45 projecting a pattern' of discrete separated point light 

more than one rectilinear sequence. sources and forming a real image of the pattern of 

3. Apparatus according to claim 1, wherein the pat- point i ight sources a t a position located in front of 
tern of discrete separated point light sources includes an the eye, 

asymmetrical shape having a plurality of lines of point selecting and collecting a reflected image of the pat- 

hght sources. 50 t ern as reflected paraxially off the cornea, and de- 

4. Apparatus according to claim 1, wherein the pat- tecting a reflected position of substantially each 
tern of discrete separated point light sources lie along an po i nt H g ht source, as reflected from the cornea, 
invertible function of distance from the optical axis of analyzing the returned, collected pattern image, in- 
the objective lens. eluding filtering noise from the collected pattern 

5. Apparatus according to claim 1, in combination 55 image using mathematical morphological transfor- 
with a surgical microscope having a standard auxiliary mations, thresholding, or fast Fourier transfonna- 
camera mount, and the ophthalmic diagnostic instru- tion techniques, 

ment being connected to the surgical microscope via comparing the filtered collected pattern image to the 

the auxiliary camera mount, with an objective lens of undistorted pattern as projected, including analyz- 

the surgical microscope serving as the objective lens of 60 ing the relative location and spatial orientation of 

the ophthalmic diagnostic instrument. the reflected point light sources as compared to the 

6. Apparatus according to claim 1, wherein the oph- pattern as projected, and 

thalmic diagnostic instrument includes means for fold- deriving mathematically a close approximation of a 

ing the pattern of discrete separated point light sources corneal surface shape that would give rise to such 

onto the optical axis of the instrument, toward the cor- 65 a collected pattern image. 

nea, with the means for projecting the pattern including 16. The method of claim 15, wherein the pattern of 

a source of the pattern off-axis from the optical axis and discrete separated point light sources comprises a gen- 

from the path of the returned, distorted pattern image. erally asterisk shaped pattern with an intersection point 



at the optical axis of the instrument, and including 
means associated with the pattern for establishing a 
readily identifiable rotational orientation of the pattern. 

17. The method of claim 15, wherein the pattern of 
discrete separated point light sources can be repre- 
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sented as any uniquely invertible function of position 
with respect to the optical axis of the image detection 
means. 

***** 
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[571 ABSTRACT 

A system, and its method, for measuring the optical prop- 
erties of an eye employs optical components for detecting 
both the reflection of a first light beam from the anterior 
surface of the eye, and the reflection of a second light beam 
from the retina of the eye. Sensors are included to receive 
and separate each of these reflected light beams into a 
plurality of individual beams, each having its own specific 
optical path length. The optical path lengths of individual 
beams reflected from the cornea are collectively used to 
create a topographical map of the cornea's anterior surface 
and the optical path lengths of individual beams reflected 
from the retina are collectively used to create an acuity map 
of the entire eye. Further, there are additional optical com- 
ponents which respectively determine the position of the 
eye, a length for the eye and an aberration for the relaxed 
lens of the eye. A computer is then used to compare the 
topographical map with the acuity map, while compensating 
for the lens aberration, to construct a topography for the 
posterior surface of the cornea. 

20 Claims, 3 Drawing Sheets 
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METHOD AND APPARATUS FOR 
MEASUREMENT OF THE REFRACTIVE 
PROPERTIES OF THE HUMAN EYE 

FIELD OF THE INVENTION 

llie present invention pertains generally to methods and 
apparatus for performing diagnostic evaluations of the 
refractive properties of a human eye. More particularly, the 
present invention pertains to the compilation of information 
and measurements that are useful for determining and select- 
ing the procedures that are necessary and appropriate for 
vision correction. The present invention is particularly, but 
not exclusively, useful for creating acuity and topographical 
maps of the refractive power of the human eye which can be 
used for the prescription of corrective elements, such as 
contact lenses and glasses, or for planning the conduct of 
refractive surgery. 

BACKGROUND OF THE INVENTION 

In the perfect eye, an incoming beam of light is focused 
through the cornea and through the crystalline lens in a way 
which causes all of the light from a point source to converge 
at the same spot on the retina of the eye. This convergence 
occurs because all of the optical path lengths, for all light in 
the beam, are equal to each other. Stated differently, in the 
perfect eye, the time for all light to transit through the eye 
will be the same regardless of the particular path that is taken 
by the light. 

Not all eyes, however, are perfect. The consequences of 
this are that light path lengths through the eye become 
distorted and are not all equal to each other. Thus, light from 
a point source that transits an imperfect eye will not neces- 
sarily be focused on the retina, or to the same spot on the 
retina. 

As light enters and passes through an eye it is refracted at 
the anterior surface of the cornea, at the posterior surface of 
the cornea, and at the surfaces of the crystalline lens. It is 
after all of these refractions have occurred that the light 
finally reaches the retina. As indicated above, in the case of 
the perfect eye, all of these refractions result in no overall 
change in the optical path lengths of light in the incoming 
beam. Therefore, any deviations which result in unequal 
changes in these optical path lengths are indicative of 
imperfections in the eye which may need to be corrected. 

In general, vision difficulties in the human eye can be 
characterized by the changes and differences in optical path 
lengths that occur as light transits through the eye. These 
difficulties are not uncommon. Indeed, nearly one half of the 
world's population suffers from imperfect visual perception. 
For example, many people are near-sighted because their 
eyeballs are "too long" (myopia). As a result, the sharp 
image of an object is generated not on the retina, but in front 
of or before the retina. Therefore, for a myopic person a 
distant scene appears to be more or less blurred. On the other 
hand, hyperopia is a condition wherein the error of refraction 
causes rays of light entering the eye parallel to the optic axis 
to be brought to a focus behind the retina. This happens 
because the eyeball is "too short" from front to back. This 
condition is commonly referred to as far-sightedness. Unlike 
the myopic person, a hyperopic, or far-sighted, person will 
see a near scene as being more or less blurred. 

Another refractive malady is astigmatism. Astigmatism, 
however, is different than either myopia or hyperopia in that 
it results from an unequal curvature of the refractive surfaces 
of the eye. With astigmatism, a ray of light is not sharply 
focused on the retina but is spread over a more or less diffuse 
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area. Further, there are even higher order refractive maladies 
of interest for vision correction which include coma and 
spherical aberration. More specifically, coma is an aberra- 
tion in a lens or lens system whereby an off-axis point object 
is imaged as a small pear-shaped blob. Coma is caused when 
the power of the zones of the lens varies with distance of the 
zone from the axis. Spherical aberration, on the other hand, 
results from loss of definition of images that are formed by 
optical systems, such as an eye. Such aberrations arise from 
the geometry of a spherical surface. 

In the past, simple refractive errors of the human eye 
(myopia, hyperopia and astigmatism) have been corrected 
conventionally with glasses, dating back to the year 1750. 
More recently, contact lenses, which were invented about 50 
years ago, have been useful for correcting these same more 
simple refractive errors. Further, refractive laser surgery 
using Excimcr UV-lascrs is receiving increased popularity. 
Thus far, however, all of these techniques for correcting 
optical impairments of the eye have been limited to the 
correction of errors from near-sightedness (myopia) or far- 
sightedness (hyperopia), and to the cylindrical refractive 
errors, the so-called astigmatism. 

As noted above, vision and its refractive errors can be 
quite complex. Similar to every other optical system, in 
addition to the simple refractive errors, the human eye also 
shows higher order refractive errors ("aberrations") such as 
coma and spherical aberration mentioned above. In all cases, 
aberrations result when an ideally flat 'wavefront' (i.e. a 
condition wherein all optical path lengths are equal) is 
distorted by a real-world optical system. In some cases, 
these distortions occur in a very complex way. In the trivial 
case, simple distortions like nearsightedness and far- 
sightedness would result in an uncomplicated bowl-like 
symmetrical distortion. With higher order aberrations, 
however, the result is a complex non-symmetrical distortion 
of the originally flat wavefront. It is these non-symmetrical 
distortions which are unique for every optical system, 
including every single person's eye, and which lead to 
blurred optical imaging of viewed scenes. 

It happens that refractive errors (aberrations or 
distortions) arc stronger when light not only passes through 
the center of an optical system, but also through the outer 
regions of the system. Specifically, these aberrations are 
more pronounced under critical lighting conditions (e.g., 
twilight). For example, it is well known that people have a 
comparably small pupil in bright daylight. As the light level 
decreases, however, the pupil becomes dilated in order to let 
more light pass through to the retina. With dilation, in 
addition to passing through the center of the eye light rays 
will also pass through the outer region of the eye (e.g. the 
optical system), where the optical quality is low. Thus, even 
persons with normal 20/20 vision have decreased visual 
acuity under critical light conditions due to increased higher 
order aberrations. 

A typical approach for improving the vision of a patient 
has been to first obtain measurements of the eye which relate 
to the topography of the anterior surface of the cornea. 
Specifically, such measurements are made to determine the 
Zernike polynomials. The Zernike polynomials are then 
used to mathematically describe and to model the anterior 
surface of the cornea. In accordance with this practice, 
depending on the order of the Zernike polynomial a certain 
refractive condition of the eye can be described. For 
example, the first order terms of the Zernike polynomials 
describe the tilt of a wavefront while second order terms 
describe myopia, hyperopia and astigmatism. Third order 
terms then describe coma and fourth order terms describe 
i.e. spherical aberration. 
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Until now, the complex aberrations of the human eye 
involving coma and spherical aberration could not be mea- 
sured and, therefore, they could not be corrected. Further, 
even today, the measurement of the 'standard' so-called 
simple refractive errors is still not fully objective. In fact, 
presently the patient's vision is usually categorized using an 
autorefractor for measuring near-sightedness, far- 
sightedness, and astigmatism. In the process, cooperation of 
the patient is crucial for obtaining even rough realistic 
results with these systems. Still, after this rough initial 
measurement, the optometrist has to use correction lenses in 
a subjective procedure to find the corrective strength that is 
best suited for the patient. To a great extent, these limitations 
have been caused by an inability to determine a topography 
for the posterior surface of the eye in addition to determining 
the topography of the anterior surface. Further, there has 
been little attention given to the peripheral areas of the 
cornea where spherical aberrations become more prominent 
as the pupil of the eye dilates. In order to overcome these 
deficiencies, it is necessary to evaluate new ways and 
methods for measuring the refractive characteristics of the 
cornea. 

Heretofore, it has been a common practice to analyze and 
describe light beams in terms of wavefronts and aberrations 
of a wave front. In this regard, the Zernike polynomials have 
been helpful. Alight beam, however, can be conceptualized 
in a different way; other than as a wavefront. It can also be 
thought of in terms of a plurality of individual beams, each 
of which has its own optical path length. Specifically, by 
way of comparison, at any particular point in time a wave- 
front can be thought of as being the temporal lengths of the 
various optical paths that have been traveled by individual 
light beams from the origin or source of the light. Thus, a 
light beam with a flat or planar wavefront is equivalent to a 
light beam wherein all light in the beam has traveled on 
optical paths that have the same temporal length. A wave- 
front can be distorted by imperfections in the eye and result 
in so-called wave aberrations. In terms of optical path 
lengths, these same aberrations can be thought of as result- 
ing from differences in the optical path lengths of individual 
beams which are caused by undesirable refractions of light 
as it passes through the eye. 

As discussed above, until now vision correction has been 
primarily concerned with reshaping the cornea using data 
that is collected about the topography of the anterior surface 
of the eye. A good example of technology that is useful for 
this purpose is provided in U.S. Pat. No. 5,062,702 which 
issued to Bille for an invention entitled "Device for Mapping 
Corneal Topography." The posterior surface of the eye, 
however, also affects the refraction of light as it passes 
through the eye. Thus, additional information about the 
thickness of the cornea is necessary for more precise refrac- 
tive corrections. To this end, a map of the posterior surface 
of the cornea would undoubtedly be useful. Further, while 
gross approximations of the lower order visual aberrations 
using the Zernike polynomials may be useful for limited 
purposes, the superficial models provided by the Zernike 
polynomials become quite cumbersome and less precise 
when higher order aberrations are concerned. 

In light of the above, it is an object of the present 
invention to provide a method and apparatus for measuring 
the refractive properties of the human eye that arc capable of 
creating a topographical map of the anterior surface of the 
eye and an acuity map of the refractive power of the entire 
human eye which are useful either for the prescription of 
corrective elements or to plan for surgery. Another object of 
the present invention is to provide a method and apparatus 
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for measuring the refractive properties of the human eye by 
considering the net effect on individual beams as they each 
pass through the eye. Yet another object of the present 
invention is to provide a method and apparatus for measur- 
ing the refractive properties of the human eye which, in 
addition to myopia, hyperopia and astigmatism, can also be 
used to determine higher order refractive error (aberrations) 
such as coma and spherical aberration. Still another object of 
the present invention is to provide a method and apparatus 
for measuring the refractive properties of the human eye 
which are effectively easy to use, relatively simple to operate 
and implement, and comparatively cost effective. 

SUMMARY OF THE PREFERRED 
EMBODIMENTS 

In accordance with the present invention, a system for 
measuring the refractive properties of the eye includes an 
optical subsystem for precisely determining the position of 
the eye. More specifically, this subsystem includes a pupil 
camera for establishing the general x-y position of the eye, 
and a confocal detector for precisely establishing the z 
position of the eye. 

Once the position of the eye has been determined and 
stabilized, a light source is activated to direct a light beam 
toward the eye for reflection from the anterior surface of the 
eye. Alenslet is also positioned in the system to separate the 
light that is reflected from the anterior surface of the eye into 
a plurality of individual beams which, depending on the 
topography of the cornea's anterior surface, will each have 
their own optical path length. These optical path lengths 
may, or may not, be equal to each other. The individual 
beams of this plurality are then directed toward a sensor, 
which cooperates with a computer to create a digital topo- 
graphical map of the cornea. This topographical map is thus 
based on the optical path lengths of the individual light 
beams that are reflected from the anterior surface of the 
cornea. For purposes of the present invention, light in the 
beam that is to be reflected from the anterior surface of the 
cornea will, preferably, have a wavelength of approximately 
840 nm and be focused to the center of curvature of the 
cornea. 

Another light source is positioned in the system to direct 
a light beam toward the eye for reflection from the retina of 
the eye. Preferably, the light in this beam is collimated as it 
arrives at the eye, and the beam has a diameter less than 
about 2 mm as it passes through the pupil toward the retina. 
The light that is then reflected from the retina fills the pupil 
and is directed toward a lenslet where, like the light reflected 
from the cornea, it is separated into a plurality of individual 
beams. Depending on the particular refractions which have 
occurred as these individual beams traverse the eye, they 
may or may not be equal to each other. The individual beams 
of this plurality are then directed toward a sensor. Like the 
sensor associated with the individual beams reflected from 
the cornea, this sensor also cooperates with the computer. 
The individual beams of light reflected from the retina, 
however, are used to create a digital acuity map of the entire 
eye. The acuity map that is thus created, is based on the 
optical path lengths of the individual light beams that are 
reflected from the retina. For purposes of the present 
invention, light in the beam that is to be reflected from the 
retina will, preferably, have a wavelength of approximately 
780 nm. 

As intended for the present invention, in order to deter- 
mine the topography of the posterior surface of the cornea, 
the refractive effects of the crystalline lens in the eye must 
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be accounted for. This can be accomplished in either of two 
ways. For the myopic eye, the patient is required to con- 
centrate on a point at infinity. The crystalline lens in the 
patient's eye will then be relaxed and, therefore, contribute 
only its basic relaxed state refraction. On the other hand, for 5 
the hyperopic eye, or the infant eye, several successive 
measurements from successive pluralities of the individual 
beams that are reflected from the retina are required. It 
happens, that these measurements, taken in defocused con- 
ditions of the lens, are collectively proportional to the 10 
relaxed condition. Thus, by using curve fitting techniques, 
the relaxed condition for the lens can be extrapolated from 
the data. A determination of the topography of the posterior 
surface of the cornea is then made by subtracting the 
topography of the anterior surface of the eye (the topo- 15 
graphical map) from the data taken for determining the 
refractive characteristics of the entire eye (the acuity map). 

BRIEF DESCRIPTION OF THE DRAWINGS 

The novel features of this invention, as well as the 20 
invention itself, both as to its structure and its operation, will 
be best understood from the accompanying drawings, taken 
in conjunction with the accompanying description, in which 
similar reference characters refer to similar parts, and in 
which: 

FIG. 1 is a schematic drawing of the system of the present 
invention; 

FIG. 2 A is a side cross sectional view of an eye shown in 
relationship with an incoming planar wavefront of light; 

FIG. 2B is a side cross sectional view of an eye shown in 
relationship with an outgoing wavefront of light that has 
been distorted by the refractive properties of the eye; 

FIG. 3 is a side cross sectional view of an eye shown with 
light reflecting from the retina; 35 

FIG. 4 is a side cross sectional view of an eye shown with 
light reflecting from the anterior surface of the cornea; and 

FIG. 5 is a graph showing the relationship between 
defocus conditions of an eye and the spherical aberration 
exhibited by the lens of the eye in a respective defocused 40 
condition. 

DESCRIPTION OF THE PREFERRED 
EMBODIMENTS 

Referring initially to FIG. 1, a system for determining the 45 
refractive properties of the human eye is shown in a sche- 
matic drawing and is generally designated 10. As impliedly 
indicated in FIG. 1, use of the system 10 is intended to 
provide diagnostic information about the eye 12. In order to 
do this, the system 10 employs four different light sources 50 
and uses four different wavelengths, all for different pur- 
poses. More specifically, the system 10 includes a light 
source 14 which is preferably a laser diode that generates a 
light beam having a wavelength of approximately eight 
hundred and forty nanometers (840 nm). Another light 55 
source 16 is provided which is preferably a laser diode that 
will generate a light beam having a wavelength of approxi- 
mately seven hundred and eighty nanometers (780 nm). 
There is still another light source 18 which is preferably a 
laser diode that will generate a light beam having a wave- 60 
length of approximately nine hundred and thirty nanometers 
(930 nm). Finally, there is an illuminator 20 which can 
include a plurality of infrared diodes that will collectively 
generate a light beam having a wavelength of approximately 
nine hundred and eighty nanometers (980). 55 

As intended for the system 10 of the present invention, a 
computer 22 is used to evaluate the light that is emitted from 



each of the above mentioned light sources 14, 16 and 18, and 
from the illuminator 20. More specifically, this evaluation is 
conducted by the computer 22 after the light from its 
respective source has been directed toward, and reflected in 
some way from the eye 12. With the objective of considering 
the collective effect of all light that is reflected in the system 
10 from the eye 12, it is perhaps best to first consider each 
source of light individually, and to track light from that 
particular source as it passes through the system 10. 

Beginning at the light source 14, as stated above, a light 
beam 24 having a wavelength of approximately eight hun- 
dred and forty nanometers (840 nm) is generated. Further, 
the light source 14 directs this beam 24 toward a dichroic 
beam splitter 26 which will pass 840 nm light, but which will 
otherwise reflect light that is substantially below 840 nm 
(e.g. 780 nm). After passing through the beam splitter 26, the 
light beam 24 is then reflected by a polarizing beam splitter 
28 for further transmission by a beam expander that includes 
the lenses 30 and 32. The light beam 24 then passes through 
a dichroic beam splitter 34 which, like the beam splitter 26, 
will pass 840 nm light but reflect 780 nm light. After passing 
through the beam splitter 34, the light beam 24 is expanded 
by the collective effect of lenses 36, 38 and 40, and it passes 
through the dichroic beam splitter 42 toward the dichroic 
beam splitter 44. For purposes of the present invention, the 
dichroic beam splitter 42 will pass light having wavelengths 
below about 900 nm and will reflect light having wave- 
lengths above about 900 nm. On the other hand, the dichroic 
beam splitter 44 will pass light having wavelengths above 
about 830 nm and will reflect light having wavelengths 
below about 830 nm. Thus, the light beam 24 will pass 
through both of the beam splitters 42 and 44. Upon passing 
through the beam splitter 44, the light in light beam 24 
passes through a quarter wave plate 46 where it is rotated 
about forty five degrees (45°). The light beam 24 is then 
focused by a moveable lens 48 onto the eye 12. 

In reflection from the eye 12, the light beam 24 passes 
back through the quarter wave plate 46 where it is again 
rotated an additional forty five degrees (45°). Thus, it is now 
rotated about ninety degrees (90°) relative to the fight in 
light beam 24 as it was being initially emitted from light 
source 14. Further, the reflected light beam 24 passes back 
through the beam splitters 44, 42, and 34. Due to its dual 
rotation by the quarter wave plate 46, however, light beam 
24 will not be reflected by the polarizing beam splitter 28. 
Instead, the polarizing beam splitter 28 will pass the light 
beam 24 toward a lenslet array 50 where the beam 24 is 
separated into a plurality of individual beams. These indi- 
vidual beams are all mutually parallel, and from the lenslet 
array 50 they are directed toward a dichroic beam splitter 52 
which, like beam splitters 26 and 34, will pass light having 
a wavelength of 840 nm. After passing through the beam 
splitter 52, the individual beams which now comprise the 
beam 24 are received by an area sensitive detector 54 and are 
then passed as a plurality of respective signals via the line 56 
to the computer 22. Preferably, the area sensitive detector 54 
is a CCD of a type well known in the pertinent art. 

In a manner similar to that set forth above for describing 
the path taken by light beam 24, consider now the light beam 
58 that is generated by the light source 16. As indicated 
above, the light beam 58 will have a wavelength that is about 
780 nm. Therefore, the light beam 58 will be reflected by the 
dichroic beam splitter 26 and the polarizing beam splitter 28. 
Unlike the light beam 24, however, the light beam 58 will be 
reflected by the beam splitter 34 and directed toward the lens 
60 and pinhole 62. A dichroic beam splitter 64 is then 
provided to direct the light beam 58 through a lens 66, 
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toward a turning mirror 68, and through a focusing lens 70. 
It can be noted that for purposes of the present invention, the 
beam splitter 64 needs to be able to pass visible light below 
a wavelength of 780 nm (i.e. light in the wavelength range 
of 780-380 nm). After being reflected by the beam splitter 5 
44, the light beam 58 is rotated by quarter wave plate 46 and 
directed by focusing mirror 48 toward the eye 12. 
Importantly, as light beam 58 is directed toward the eye 12, 
the light in light beam 58 will be substantially collimated 
and have a beam diameter that is approximately two milli- 
meters (2 mm). 

When reflected from the eye 12, the light beam 58 will 
again be rotated by the quarter wave plate 46. Thus, it will 
be reflected by the beam splitter 44, turned by turning mirror 
68 and reflected by the beam splitters 64 and 34. Like the 15 
light beam 24, light beam 58 will be passed by the polarizing 
beam splitter 28 toward the lenslet array 50. Also like the 
light beam 24, the light beam 58 will be separated into a 
plurality of individual light beams by the lenslet array 50. 
The plurality of individual light beams which now comprise 9Q 
the light beam 58 are reflected by the beam splitter 52 and 
are directed toward an area sensitive detector 72 where the 
individual beams arc converted into respective signals for 
further transmission via line 74 to the computer 22. 

The light source 18, as mentioned above, will generate a 2 5 
light beam 75 which has a wavelength of approximately 930 
nm. As shown in FIG. 1, this light beam 75 will pass through 
a polarizing beam splitter 76, a lens 78 and a dichroic beam 
splitter 80, before being reflected toward the eye 12 by the 
beam splitter 42. Importantly, as the beam 75 is directed 30 
toward the eye 12 it will pass through and be rotated by the 
quarter wave plate 46. The light beam 75 is then reflected 
from the eye 12, and light in the reflected light beam 75 will 
again pass through, and be rotated by, the quarter wave plate 
46. At the beam splitter 42, the light beam 75 will be directed 35 
back toward the polarizing beam splitter 76. This time, 
however, light beam 75 does not pass through the beam 
splitter 76. Instead, due to its rotations by the quarter wave 
plate 46, the light beam 75 is reflected by the polarizing 
beam splitter 76 through a pinhole 81 toward the confocal 40 
detector 82. A signal that is generated by the confocal 
detector 82 in response to its reception of the light beam 75 
is then passed via line 84 to the computer 22. 

As indicated above, the illuminator 20 generates a light 
beam 86 which has a wavelength of about 980 nm. For the 45 
present invention, the illuminator 20 can include either a 
plurality of separate infrared diodes, or it can be configured 
as a ring. In either case, as shown in FIG. 1, the resultant 
light beam 86 is pointed directly at the eye 12. Upon 
reflection of the light beam 86 from the eye 12, FIG. 1 50 
indicates that the beam passes through the beam splitter 44 
but is reflected by both the beam splitter 42 and the beam 
splitter 80. Specifically, insofar as the beam splitter 80 is 
concerned, it will reflect light such as the light in light beam 
86 which has a wavelength greater than about 950 nm. 55 
Accordingly, the light in light beam 86 which has been 
reflected from the eye 12 will be received by a pupil camera 
88 and a responsive signal generated by the pupil camera 88 
will be sent via line 90 to the computer 22. 

FIG. 1 also shows that the computer 22 is connected via 60 
a line 92 with the lens 70. With this connection the computer 
22 is able to adjust the focus of light beam 58 that is 
provided by lens 70. Further, FIG. 1 shows that the computer 
22 is connected via a line 94 with the lens 48. With this 
connection the computer 22 is able to adjust the focus of 65 
light beams 24 and 75. Also, FIG. 1 shows that the computer 
22 can include a frame grabber 96 which will provide visual 



displays of the signals that are received from the area 
sensitive detectors 54 and 72, as well as the signals that are 
received from the confocal detector 82 and the pupil camera 
88. 

OPERATION 

In overview, for the purposes of the operation of system 
10, light can be conceptualized in either of two ways. Firstly, 
light can be thought of in terms of wavefronts. Secondly, 
light can be thought of as being a collective bundle of many 
different individually separate beams. These two concepts, 
of course, must be related to each other if they are to 
describe the same light beam. Accordingly, in order to 
reconcile one concept with the other, a wavefront can be 
thought of as being a spatial representation of the optical 
path lengths that have been traveled from a common origin 
(light source) by all of the different individual beams, at any 
given point in time. Thus, it is the case with unrefracted light 
which has traveled from a light source in the direction of 
arrow 98, as shown in FIG. 2A, that the light will exhibit a 
planar wavefront 100. Stated differently, the optical path 
length of an individual beam in this light that has traveled 
from the source to a position 102 in the wavefront 100 will 
have the same length as the optical path length of an 
individual beam that has traveled from the source to a 
position 104 in the wavefront 100. As the light in wavefront 
100 passes through the eye 12, however, the individual light 
beams will be refracted differently. 

Anatomically, it is necessary for light entering the eye 12 
to pass through the cornea 106 and the crystalline lens 108 
of the eye 12 before coming into contact with the retina 110. 
It is known that our visual perception is dependent on how 
this light comes into contact with the retina 110 and, it is 
known, of course, that in accordance with Snell's law this 
light will be refracted by the cornea 106 and the lens 108 as 
it passes into the eye 12. Further, whatever light there is that 
is reflected from the retina 110 to pass back through the eye 
12 and away from the eye 12 will also be refracted by the 
lens 108 and the cornea 106. The result of all this refraction 
may likely be a distorted wavefront 112 which is traveling 
away from the eye 12 in the direction of arrow 114, such as 
shown in FIG. 2B. By comparing FIG. 2Awith FIG. 2B it 
will be noted that, due to refractions caused by the eye 12, 
the optical path length of the individual beam which has 
traveled from position 102 in the planar wavefront 100 to the 
position 102' in the distorted wavefront 112 will be different 
from the optical path length of the individual beam which 
has traveled from position 104 to position 104'. As appre- 
ciated by the present invention, the differences in these 
optical path lengths is indicative of the respective refractions 
that were experienced by the individual beams as they 
transited the eye 12. 

In the operation of the system 10 of the present invention, 
it may first be necessary to calibrate the system 10. This can 
be done by replacing the eye 12 with a flat mirror (not 
shown). Light can then be sequentially passed through the 
system 10 from the light sources 14, 16 and 18 and from the 
illuminator 20 for reflection back through the system 10 
from the flat mirror. In this way, signals can be generated 
from the area sensitive detectors 54 and 72, from the 
confocal detector 82 and from the pupil camera 88. The 
signals which are thus generated will be indicative of 
inherent optical aberrations in the system 10 and can sub- 
sequently be used to compensate actual signals generated by 
the eye 12. 

Once the system 10 has been calibrated, it is desirable to 
determine an exact spatial position for the eye 12 in x-y-z. 
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This is done by using the confocal detector 82 and the pupil 
camera 88. Specifically, in order to establish a "z" position 
for the eye 12, the light source 18 is activated to generate the 
light beam 75. For the present invention, the light beam 75 
is focused by the lens 48 to obtain a specular reflection of the 
light beam 75 from the apex 114 of the cornea 106 (see FIG. 
3). Depending on the position of the lens 48, which is sensed 
by the computer 22, when this specular reflection is obtained 
the position of the eye 12, and more specifically the position 
of the apex 114 of eye 12, in a "z" direction along the visual 
axis of the eye 12 is established. In order to establish an 
"x-y" position for the eye 12, the illuminator 20 needs to be 
activated. Specifically, using intensity differences in the 
reflection of light beam 86 from the eye 12, as sensed by the 
pupil camera 88 and ultimately by the computer 22, the 
"x-y" position of the eye 12 is established. For the present 
invention, the intensity differences used for this measure- 
ment are cause by the contrast between the iris 116 and the 
lens 108 at the periphery of the pupil 118. While the "z" 
position has been considered as having been taken first in 
this discussion, it will be appreciated by the skilled artisan 
that the "x-y" determination may, in fact, be made first. 

Once the system 10 is calibrated and the position of the 
eye 12 in "x-y" and "z" has been established, refractive 
measurements of the eye can be made. In light of the above 
discussion, it will be noted that these measurements need to 
be considered with knowledge of the length of the eye 12, 
and with an understanding that the crystalline lens 108 is in 
its basic relaxed state of refraction. Insofar as a measurement 
for the length of the eye 12 is concerned, this can be done 
by activating the light source 114 while the eye 12 is focused 
to a point 124 at infinity (see FIG. 1). The focusing lens 70 
is then moved as required by the computer 22, or manually 
by the operator of system 10, until the light beam 24 from 
light source 114 is focused onto the retina 110. By using the 
position of the lens 70 for this focal condition, the computer 
22 is able to establish a position of the retina 110. Then, by 
knowing the position of the retina 110, and by knowing the 
location of the apex 114 that was obtained from earlier 
measurements of the "z" position of the eye 12, the length 
of the eye 12 can be determined. This measurement, of 
course, will also determine whether the eye 12 is myopic or 
hyperopic. 

Referring to FIG. 3, it will be appreciated that if the light 
beam 24 is focused by the eye 12, without any correction, to 
a location 120 in front of the retina 110, and movement of 
the lens 70 is necessary to move the focal point of light beam 
24 from the location 120 backward toward the retina 110, the 
eye 12 is myopic. On the other hand, if the light beam 24 is 
focused by the eye to a location 122 behind the retina 110, 
again without any correction, and movement of the lens 70 
is necessary to move the focal point of light beam 24 from 
the location 122 forward toward the retina 110, the eye 12 
is hyperopic. The determination of whether the eye 12 is 
myopic or hyperopic is important, not only in its own right, 
but also for subsequent refractive measurements. 
Importantly, while the myopic eye 12 focuses on the point 
124 at infinity, it can be taken that the crystalline lens 70 will 
be in its basic relaxed state of refraction. On the other hand, 
as mentioned above and as explained in more detail below, 
when the eye 12 is hyperopic of is the eye of an infant 
several successive measurements need to be taken. 

The initial measurement for the general topography of the 
cornea 106 of the eye 12 is made using the light beam 24 that 
is generated by light source 14, and which has a wavelength 
of 840 nm. By cross referencing FIG. 1 with FIG. 4 it will 
be seen that the light beam 24 is focused by the lens 48 
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toward the center of curvature 126 of the cornea 106. In 
doing this, the light beam 24 covers a distance 128 on the 
anterior surface 130 of the cornea 106 that is equal to about 
seven millimeters (7 mm). Importantly, some of the light in 
the light beam 24 will be reflected by the anterior surface 
130 of the cornea 106 and will be directed back through the 
system 10 to the lenslet 50 where it is separated into a 
plurality of individual beams. These individual beams in the 
reflected light of light beam 24 are then detected by the area 
sensitive detector 54 which generates signals that are sent 
via line 56 to the computer 22. Using these particular 
signals, the computer 22 is able to create a topography map 
of the anterior surface 130 of the cornea 106. 

For the present invention, refractive measurements of the 
entire eye 12 are made using the light beam 58 that is 
generated by light source 16. The shorter wavelength of 780 
nm is selected for the light beam 58 because it is near the 
visible range and it will, therefore, more easily travel 
through the eye 12 than will the longer wavelength light 
beams used for other purposes in the system 10. It is an 
important consideration that the light beam 58 have a 
relatively small cross section as it initially enters the eye 12. 
This is so in order to minimize light refractions that arc 
caused as the light beam 58 travels through the eye 12 
toward the retina 110. For the present invention, the light 
beam 58 is, preferably, confined to about two millimeters (2 
mm) diameter. Also, the light beam 58 is adjusted by the 
optics along its beam path so that as the light beam 58 leaves 
the lens 48 and travels toward the eye 12, it is substantially 
collimated when it arrives at the anterior surface 130 of the 
cornea 106. 

Returning to FIG. 3, it will be seen that the light in light 
beam 58 will be reflected from the retina 110 as a light beam 
58' which fills the pupil 118. This reflected light beam 58' 
then is passed back through the system 10 to the lenslet 50 
where, like the beam 24, it is separated into a plurality of 
individual beams. Also like the individual beams of light 
beam 24, the individual beams of light beam 58 are passed 
to an area sensitive detector 72 where signals are generated 
for transmission via line 74 to the computer 22. More 
specifically, the individual beams of light beam 58 are 
collectively used to generate an acuity map of the eye 12 
which is indicative of the refractions experienced by light as 
it passes through the eye 12. 

For a myopic eye 12, which will remain in its basic 
relaxed state of refraction while focused to infinity on the 
point 124 (see FIG. 1), a topography for the posterior surface 
132 of the cornea 106 can be determined using computer 22. 
Basically, this is done by subtracting the topography map 
data for the anterior surface 130 of the eye 12, and the basic 
relaxed state refraction for the crystalline lens 108, from the 
acuity map of the entire eye 12. The result is data which can 
be used directly to determine the topography for the poste- 
rior surface 132 of the eye 12. 

For a hyperopic eye 12 or an infant eye 12, which is not 
able to establish its basic relaxed state of refraction while 
focused to infinity on the point 124, successive measure- 
ments need to be taken and the collected data extrapolated 
to determine the basic relaxed state of refraction. This is 
possible because it is known that there is a generally linear 
relationship between each defocus condition of the crystal- 
line lens 108 in the eye 12, and the corresponding spherical 
aberrations caused by the lens 108 (sec FIG. 5). Therefore, 
by taking a series of successive measurements for the acuity 
map (i.e. using the light beam 58 from light sources 16) a 
plurality of data points 134 (of which the data points 134a, 
1346 and 134c are representative), can be plotted. In FIG. 5 
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it is seen that the plots of data points 134 can be used to 
identify a line 136 and that the point 138 can be extrapolated 
and be considered equivalent to the conditions extant when 
the crystalline lens 108 is in its basic relaxed state of 
refraction. In a manner as disclosed above in consideration 
of the myopic eye 12, this data can be used with the 
topography map of the anterior surface 130 of the cornea 
106 and the acuity map of the entire eye 12 to determine a 
topography map for the posterior surface 132 of the eye 12. 
In any case, all of the data collected will give the operating 
physician a much more detailed measurement of the 
anatomy of the eye 12 which will be useful for the prescrip- 
tion of corrective elements or for planning the conduct of 
refractive surgery. 

While the particular method and apparatus for measure- 
ment of the refractive properties of the human eye as herein 
shown and disclosed in detail is fully capable of obtaining 
the objects and providing the advantages herein before 
stated, it is to be understood that it is merely illustrative of 
the presently preferred embodiments of the invention and 
that no limitations are intended to the details of construction 
or design herein shown other than as described in the 
appended claims. 

What is claimed is: 

1. A system for measuring optical properties of an eye, the 
eye having, in order, a cornea, a lens and a retina, said 
system comprising: 

an optical means for directing a first light beam for 
reflection from the cornea of the eye as a reflected first 
light beam; 

a lenslet array for separating said reflected first light beam 
into a plurality of first individual beams with each said 
first individual beam having an optical path length; 

an optical means for aiming a second light beam through 
the cornea, and through the lens, for reflection from the 
retina of the eye as a reflected second light beam; 

a lenslet array for separating said reflected second light 
beam into a plurality of second individual beams with 
each said second individual beam having an optical 
path length; 

a computer means for using said plurality of first indi- 
vidual beams to create a topographical map of the 
cornea indicative of the optical path lengths of said first 
individual beams, and for using said plurality of second 
individual beams to create an acuity map of the eye 
indicative of the optical path lengths of said second 
individual beams; and 

a comparator means operable with said computer means 
for comparing said topographical map of the cornea 
with said acuity map of the eye to determine selected 
optical properties of the eye. 

2. A system as recited in claim 1 wherein the optical 
properties of the eye are determined by differences between 
the optical path lengths of said first individual beams used 
for creation of said topographical map, and by differences 
between the optical path lengths of said second individual 
beams used for creation of said acuity map. 

3. A system as recited in claim 1 further comprising a lens 
for focusing said second light beam onto the retina to 
determine a length for the eye. 

4. A system as recited in claim 1 wherein said computer 
means analyzes said topographical map to model an anterior 
surf ace of the cornea. 

5. A system as recited in claim 1 wherein said topographi- 
cal map is indicative of aberrations for the anterior surface 
of the cornea, wherein said acuity map is indicative of 
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aberrations for the entire eye, and wherein said comparator 
means separates the aberrations of the anterior surface of the 
cornea and a predetermined aberration for the lens, from the 
aberrations of the entire eye to determine aberrations for the 
posterior surface of the cornea. 

6. Asystem as recited in claim 5 wherein said aberrations 
of the posterior surface of the cornea are used by said 
computer means to model the posterior surface of the 
cornea. 

7. A system as recited in claim 5 wherein the spherical 
aberration of the lens is determined with the lens in a relaxed 
state. 

8. A system as recited in claim 5 wherein the predeter- 
mined aberration of the lens is determined by said computer 
means using a plurality of sequentially created said acuity 
maps. 

9. A system as recited in claim 1 further comprising: 

an optical means for sending a third light beam for 
reflection from the apex of the cornea as a reflected 
third light beam; 

an illuminator for generating a pupil reflection from the 
eye; and 

a detector means for analyzing said pupil reflection for 
determining an "x-y" position for the eye, and for 
analyzing said reflected third light beam for establish- 
ing a "z" position for the eye. 

10. A system as recited in claim 1 wherein the cornea has 
a center of curvature and said first light beam is focused onto 
the center of curvature of the cornea, and wherein said first 
light beam has a wavelength equal to approximately eight 
hundred and forty nanometers (840 nm). 

11. A system as recited in claim 1 wherein the eye has a 
pupil and said second light beam is aimed toward the eye as 
collimated light for subsequent focus by the cornea and the 
lens of the eye toward the retina, and wherein said second 
light beam has a wavelength equal to approximately seven 
hundred and eighty nanometers (780 nm), and further 
wherein said second light beam has a diameter of approxi- 
mately 2 mm when entering through the pupil and a diameter 
of approximately 6 mm when emerging from the eye 
through the pupil after reflection from the retina. 

12. A system for measuring optical properties of the eye, 
the eye including a cornea having an anterior surface and a 
posterior surface, a lens and a retina, said system compris- 
ing: 

means for creating a topographical map of the anterior 
surface of the eye, said topographical map including a 
plurality of individual refractive measurements; 

means for generating an acuity map for the entire eye, said 
acuity map comprising a plurality of individual refrac- 
tive measurements; 

means for determining an aberration for a relaxed lens; 
and 

computer means for comparing said topographical map 
with said acuity map while compensating for the aber- 
ration of the lens to construct a topography for the 
posterior surface of the eye. 

13. Asystem as recited in claim 12 wherein said means for 
creating a topographical map comprises: 

an optical means for directing a first light beam for 
reflection from the anterior surface of the cornea of the 
eye as a reflected first light beam; and 

a lenslet array for separating said reflected first light beam 
into a plurality of first individual beams with each said 
first individual beam having an optical path length and 
each said optical path length being indicative of a 
respective said refractive measurement. 
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14. A system as recited in claim 13 wherein said means for 
generating an acuity map comprises: 

an optical means for aiming a second light beam through 
the cornea, and through the lens, for reflection from the 
retina of the eye as a reflected second light beam; 

a lenslet array for separating said reflected second light 
beam into a plurality of second individual beams with 
each said second individual beam having an optical 
path length and each said optical path length being 
indicative of a respective said refractive measurement. 

15. A system as recited in claim 14 further comprising: 
an optical means for sending a third light beam for 

reflection from the apex of the cornea as a reflected 
third light beam; 
an illuminator for generating a pupil reflection from the 
eye; and 

a detector means for analyzing said pupil reflection for 
determining an "x-y" position for the eye, and for 
analyzing said reflected third light beam for establish- 
ing a "z" position for the eye. 

16. A system as recited in claim 14 wherein the cornea has 
a center of curvature and said first light beam is focused onto 
the center of curvature of the cornea, and wherein said first 
light beam has a wavelength equal to approximately eight 
hundred and forty nanometers (840 nm), and further wherein 
the eye has a pupil and said second light beam is aimed 
toward the eye as collimated light for subsequent focus by 
the cornea and the lens of the eye toward the retina, and 
wherein said second light beam has a wavelength equal to 
approximately seven hundred and eighty nanometers (780 
nm) and has a diameter of approximately 2 mm when 
entering through the pupil and a diameter of approximately 
6 mm when emerging from the eye through the pupil after 
reflection from the retina. 

17. A method for measuring optical properties of the eye, 
the eye including a cornea having an anterior surface and a 
posterior surface, a lens and a retina, said method compris- 
ing the steps of: 

creating a topographical map of the anterior surface of the 
eye, said topographical map including a plurality of 
individual refractive measurements; 
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generating an acuity map for the entire eye, said acuity 
map comprising a plurality of individual refractive 
measurements; 
determining an aberration for a relaxed lens; and 
comparing said topographical map with said acuity map 
while compensating for the aberration of the lens to 
construct a topography for the posterior surface of the 
eye. 

18. A method as recited in claim 17 wherein said step of 
creating a topographical map comprises the steps of: 

directing a first light beam for reflection from the anterior 
surface of the cornea of the eye as a reflected first light 
beam; and 

separating said reflected first light beam into a plurality of 
first individual beams with each said first individual 
beam having an optical path length and each said 
optical path length being indicative of a respective said 
refractive measurement. 

19. A method as recited in claim 18 wherein said step of 
generating an acuity map comprises the steps of: 

aiming a second light beam through the cornea, and 
through the lens, for reflection from the retina of the 
eye as a reflected second light beam; and 

separating said reflected second light beam into a plurality 
of second individual beams with each said second 
individual beam having an optical path length and each 
said optical path length being indicative of a respective 
said refractive measurement. 

20. A method as recited in claim 17 further comprising the 
steps of: 

sending a third light beam for reflection from the apex of 
the cornea as a reflected third light beam; 

illuminating the eye to generate a pupil reflection from the 
eye; and 

analyzing said pupil reflection for determining an "x-y" 

position for the eye; and 
analyzing said reflected third light beam for establishing 

a "z" position for the eye. 
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Photon noise and atmospheric noise in active optical systems 
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A general theoretical analysis is made of the performance of optical systems that are designed to 
give diffraction-limited images of astronomical objects by compensating effects of atmospheric seeing 
in real time. The heart of any such system is a feedback algorithm, which expresses the controlled 
displacements of mirror surfaces as functions of the output of optical sensors. The statistical behavior 
of the system is calculated, assuming the feedback algorithm to be linear but otherwise unrestricted. 
The effect of feedback in diminishing atmospheric noise while amplifying photon noise is worked but 
in detail, taking into account photon-photon correlations. A figure of merit for system performance, 
namely a statistical average of a positive quadratic function of phase errors over mirror surfaces, is 
arbitrarily chosen. By use of this figure of merit, the optimum feedback algorithm for any given 
optical system is explicitly determined.' The optimum algorithm is independent of the quadratic 
function of phase errors that is chosen as the figure of merit. Applications of the general theory to 
particular optical systems are briefly discussed. In principle, systems optimized in this way should be 
able to give images of arbitrarily high resolution of astronomical objects brighter than about 
magnitude 14. 



The phrase active optical system in this paper means a 
system designed to provide high- resolution images of 
celestial objects, with movable optical surfaces that 
compensate in real time the effects of fluctuating atmo- 
spheric perturbation of the incoming light. Such a sys- 
tem was first proposed by Babcock, 1 and has been re- 
cently investigated experimentally by Muller and Buf- 
fington, 2 Bridges et al* , 3 Hardy et al. , 4 Miller et al. , 5 
Ogrodnik, 6 and Pearson et al. 7 Every active optical 
system consists conceptually of six parts: A, one or 
more light- collecting mirrors, either in the form of a 
conventional telescope, an interferometer, or an array 
of interferometers; B, a secondary flexible mirror, in 
which each small patch of surface can be moved inde- 
pendently to introduce a controlled variation of optical 
path length of the light reflected from that patch; C, a 
set of electromechanical relays that control the move- 
ments of the patches of surface in B; D, a television 
monitor or array of photon detectors that record the 
places and times of arrival of photons that pass through 
the system; E, a computer that accepts input from D and 
which, on the basis of this information, feeds output to 
C, so that the movement of the mirror surface B is re- 
sponsive to the optical performance of the system at 
previous times; F, a computer program or algorithm 
that instructs the computer E how to process the input 
and compute the output. The essential task that the com- 
puter program has to perform is to tell the flexible mir- 
ror what to do in order to improve the optical quality of 
the image, using for this purpose the information pro- 
vided by the image itself. 

The five components A-E of the system are hardware 
items, which have been extensively discussed by the ex- 
perimenters. 2 " 7 In this paper, we take the hardware 
components as given and concentrate our attention on the 
design of the program F. To design an effective com- 
puter algorithm is equivalent to solving a prediction 
problem, namely to predict the behavior of atmospheric 
perturbations at time t when we know only the optical 
response of the system to perturbations that occur at 



times earlier than t. The possibility of an effective pre- 
diction algorithm is, in principle, limited by two sources 
of noise, first the statistical fluctuations of the atmo- 
spheric perturbations (atmospheric noise), and second 
the statistical fluctuations of the output of the optical 
sensors (photon noise). The feedback loop, from optical 
sensor to computer to servomechanism to flexible mir- 
ror to optical sensor, has the primary function of di- 
minishing the atmospheric noise, but the photon noise 
is necessarily amplified by the feedback as the atmo- 
spheric noise is attenuated. It is to be expected that an 
optimum strength of feedback will be obtained when the 
two sources of noise are made approximately equal. 
The exact optimization of this tradeoff between atmo- 
spheric and photon noise will depend in a sensitive way 
upon the frequency spectrum of the atmospheric pertur- 
bations. 

This paper consists of three parts. In Sec. I, we 
make no restriction on the nature of the feedback algo- 
rithm except that it is to be linear and causal. That is 
to say, we let the computer outputs (mirror- patch dis- 
placements at time t) be arbitrary linear functions of the 
computer inputs (optical sensor responses) at times pre- 
vious to t. We then compute the statistical behavior of 
the whole system responding to atmospheric perturba- 
tions of given statistical characteristics. The crucial 
feature of this analysis is an exact treatment of photon- 
photon correlations. Photon-photon correlations arise 
because the probability of detection of a photon in a sen- 
sor at time t 2 is significantly changed by the mirror dis- 
placement induced by the feedback response to a single 
photon detected at an earlier time t x . The degradation 
of system performance by photon noise depends in an 
essential way on these photon- photon correlations. The 
end result of Sec. I is a compact formula (42) that ex- 
hibits the dependence of the system performance upon 
the atmospheric-fluctuation spectrum, the optical char- 
acteristics of the system, and the feedback algorithm. 

In Sec. II, we use the results of Sec. I to solve exactly 
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the problem of optimizing the feedback algorithm for a 
system whose hardware components are given. Here, 
optimization has to be defined in some arbitrary fashion. 
We choose to define optimization as the minimization of 
a mean-square deviation (averaged over both atmospher- 
ic and photon fluctuations) of optical phase errors over 
the mirror surfaces. Fortunately, it turns out that the 
optimum feedback algorithm so defined is independent 
of which particular mean-square average deviation we 
choose to minimize. The optimized algorithm may 
therefore be expected to be useful even under conditions 
in which our arbitrarily chosen criterion of optimization 
is inappropriate. 

The characteristics of the optimum algorithm are de- 
termined in two steps. First, we determine the space 
dependence of the feedback by means of a duality theo- 
rem (Theorem 1), which states that the feedback from 
optical sensor j to mirror patch k should be proportional 
to the optical response of the sensor j to a displacement 
of the patch k. Second, we determine the time depen- 
dence of the feedback by a nonlinear integral equation 
(Theorem 2), which relates the feedback directly to the 
space-time correlation structure of the atmospheric 
fluctuations. The two theorems together determine the 
optimum algorithm uniquely. It turns out that the non- 
linear integral equation (85) is identical to the Gelfand- 
Levitan equation, 8 which enables an unknown scattering 
potential to be reconstructed from the scattering phase 
shifts that it induces upon waves passing through it. It 
is perhaps not accidental that the same equation arises 
in the present context, because we are dealing with the 
reconstruction of an unknown atmospheric- scattering 
medium from observations of phase shifts in the light 
that comes through it. Nevertheless, the analogy is far 
from exact. Gelfand and Levitan were not discussing 
an optimization problem, and their space coordinate 
plays the role that in our analysis is played by the time 
coordinate. The fact that precisely the same equation 
appears in these two different physical contexts indicates 
a connection whose exact logical basis remains to be 
explored. 

In Sec. in, we indicate how the general results of 
Sees. I and II are to be applied to particular optical sys- 
tems, and make a rough estimate of the brightness of 
objects for which effective compensation of atmospheric 
perturbations might be achieved. Section III is purpose- 
ly made brief, because an informed discussion of the 
practical aspects of active optical systems lies beyond 
the competence of the author. 

I. PHOTON-PHOTON CORRELATIONS 

First, we introduce notations designed to characterize 
an active optical system in a general abstract fashion. 
The advantage of these abstract notations is that the 
equations become so compact that their structure is 
easily grasped. 

The patches of movable mirror surface are indexed 
by the letter k. It is convenient to think of k as a dis- 
crete variable, but the theory applies unchanged with a 
continuously variable k if the movable mirror is regard- 
ed as a continuously deformable surface. The photon- 



detector channels are indexed by the letter j. Here j 
must in reality be a discrete variable, but it may be con- 
venient to idealize the detector system by imagining j to 
indicate the precise position of a photon arriving at the 
image plane, or the precise wave number of a photon 
passing through a spectrometer, in which case j would 
take a continuous range of values. When we refer to the 
state of mirror patch k at the time t ki we use the letter 
y instead of the pair (k, t k ). Thus, the imposed displace- 
ment of the mirror surface is defined by a function 



(l) 



whose values are controlled by the computer output. 
Similarly, the letter x denotes a detector channel j, to- 
gether with a time t j9 so that 

P=P(x) = P(j,t s ) ( 2 ) 

denotes the probability per unit time that a photon will 
arrive at channel ; at time t $ . Because of the quantum 
nature of light, the actual photon count in channel j is 
not P but 



I = l(x)=I(j, tj). 



We have 



(3) 



(4) 



where the bracket () p denotes an average over photon 
statistics, and P is the irradiance calculated from clas- 
sical electromagnetic theory, given any particular con- 
figuration of the atmosphere and the mirror displace- 
ments. 

We suppose that the unknown atmospheric perturbation 
whose effect the system is designed to compensate is de- 
noted by a function 



a = a(v) = a(^,4) J 



(5) 



representing the fluctuation of optical path length of the 
pencil of light rays that arrives from the object at the 
mirror patch k. The field of view of the system is as- 
sumed to be so small that all rays from the object to a 
given mirror patch pass through the same isoplanatic 
patch in the sky and are subject to the same atmospheric 
perturbation. The path-length perturbation a is supposed 
to be compensated in part by the mirror displacement 
\b, so that there is a net path-length error 



s-a+b~s(k,t k ) 



(6) 



in the light waves reflected from mirror patch k. The 
statistical characteristics of the atmosphere are defined 
by a space- time correlation function 



U=U{y l ,y 2 ) = U{k l , fe 2 , h-t z ) 
= (a(k lt tx)a(k 2 , tz)) a , 



(7) 



where the bracket <> fl denotes an average over atmo- 
spheric statistics. The numerical values of U will de- 
pend on seeing conditions, wind speeds and other mete- 
orological variables that are not explicitly represented 
inEq, (7). 

When a given astronomical object is observed by the 
system, classical electromagnetic theory predicts a 
distribution of irradiance P in the various detector chan- 
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nels, with P depending in a calculable way on the path- 
length errors s given by Eq. (6). We assume that the 
feedback system is successful in keeping the errors 5 
small compared with an optical wavelength, so that the 
dependence of P on s is approximately linear. We have 
then a linear equation, - 

PU*)=Po(i) + £*i*s(M), (8) 

k 

which we write in operator notation as 

P=P Q + Bs. ^ 

Here P Q is the distribution of irradiance in the detector 
channels in the absence of atmospheric perturbations, 
and B is the optical- response kernel 



(10) 



B(x,y)=B jk 5(t j 

that gives the instantaneous optical response of the de- 
tector j to a displacement of the mirror patch k. Both 
p 0 and B depend on only the optical characteristics of 
the object and the observing apparatus. They are inde- 
pendent of the atmosphere and of the feedback system. 

The feedback algorithm is assumed to be a linear re- 
lation 



Hk, t) = £ £ dt'N kS {t- t')I{j, t'), 



or in operator notation 
6 = jN7, 



(H) 



(12) 



<s>, = [i-ira]- l fl, < 16 ) 

(D^Pv + ll-BNY'Ba. (") 

The operators on the right-hand side of Eqs. (16) and 
(17) can be inverted in a straightforward manner, by 
using Fourier transforms or otherwise. The equations 
show that the feedback reduces the mean path-length er- 
rors (s) p to a small fraction of the atmospheric pertur- 
bations a, provided that the product NB is made large 
and negative. Only the effects of photon noise set an up- 
per limit to the strength of negative feedback that can be 
usefully applied. Photon-noise effects appear in the 
two-point functions or averages of quantities bilinear in 
the photon counts. We now proceed to calculate these. 



where the feedback kernel N relates the mirror displace- 
ment b to the photon count / at earlier times. It is im- 
portant that rather than P, appears on the right-hand 
side of Eq. (12), so that photon fluctuations are included 
in the feedback. The system should be in equilibrium 
with zero feedback when the irradiance has the unper- 
turbed distribution P 0 , so that we require the feedback 
kernel to satisfy the condition 

JVP o =0. ( 13 > 

It will turn out that Eq. (13) is satisfied automatically 
by the optimized feedback kernel, which we shall com- 
pute in Sec. H. 

We have now a complete abstract definition of the sys- 
tem in terms of optical- sensor output/, mirror dis- 
placement |&, atmospheric path-length perturbation a, 
atmospheric correlation function U, unperturbed irradi- 
ance P 0 , optical response kernel B, and feedback kernel 
N. These quantities are related by the closed system of 
Eqs. (4), (6), (7), (9), (12), and (13). It remains to cal- 
culate the performance of the system by solving these 
equations. 

We begin by calculating the one-point functions or av- 
erages of quantities linear in the photon counts. These 
quantities obey linear equations obtained by averaging 
Eqs. (9) and (12) over the photon statistics, 

(s) p =a+N{(I) p -P Q ), 

(I) p -P 0 = B(s) p . 
These equations have the solution 



(14) 
(15) 



We use the abbreviated notation 



(18) 



to denote the value of the quantities s and J at places and 
times labelled by the index 1. The two-point functions 
are then 



(19) 



and similarly (s^p and (I x l£ p . Two equations for these 
functions are obtained immediately by multiplying Eq. 
(12) by the value of s or I at a second place and taking 
the expectation value. In this way, we obtain 



(20) 
(21) 



Here the notation N x means the kernel N operating on the 
space-time coordinates labelled by the index 1 only, and 
similarly for N z . Equations (20) and (21) hold, irrespec- 
tive of the time sequence of the places 1 and 2. The 
multiplier function [s x in Eq. (20), I z in Eq. (21)] plays 
only the role of a spectator and does not disturb the op- 
eration of the feedback kernel. 

The third equation for the two-point functions must be 
derived from Eq. (9). The situation is here more com- 
plicated, because the quantity on the left-hand side of 
Eq. (9) is the expectation value (I) p and not J itself. K 
we multiply Eq. (9) by the value of / at a second place, 
we obtain 

(PJJ P = P Q i(tip + *i(sih)p. ( 22) 
Now, according to Eq. (4) 

<Px/ 2 >* = «/M>*> (23) 
where the inner average is over photon statistics at the 
time fx only, whereas the outer average is over photon 
statistics at all times. It is permissible to write 

Wz)p-(h^p (24) 
if and only if t x >t 2 , so that the photons detected at time 
t x cannot have any causal influence upon photons detected 
at t 2 . Equation (22) thus becomes 

valid for t t > t z only. The function (I X I Z ) P must also con- 
tain a singular contribution that arises when the places 
1 and 2 coincide in space and time. To evaluate this 
singular contribution, we consider the quantity 
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n= [ /(. 



x) dx, 



(26) 



where the x integration is over any subset of detector 
channels and over a very short interval of time St. If 
6/ is short enough, at most one photon will be detected 
within the integration volume, and so n is a random vari- 
able taking only the values 0 and 1. We therefore have 

(n% = (n) p , (27) 

or explicitly 

f f d Xl dxz(hQ p = f dx(I) p , (28) 

to first order in 6t. Equation (28) implies that (/ 1 / 2 ) 
contains the singular contribution 

A 12 = &uP u (29) 

where 

** = QUuh)*ih-h) ( 3 °) 
is a delta function that identifies the places 1 and 2 in 
both channel index and time. Equation (25), supple- 
mented by the singular term, becomes 

(hh)p = Poi(h), + *i(sil z ) p +A 12 , (31) 
valid now for ^S:^. 



It is convenient to introduce the reduced two -point 
functions defined by 



( s l S z) Pr - ( s l s z)p - ( s l)p( s z)p, 



(32) 



and similarly for (s 1 / 2 )/ >r and ( I ih)pr- Tne reduced func- 
tions represent purely the effects of photon-photon cor- 
relations. Equations (20), (21), and (31) give three 
equations for the reduced functions, 



(SiS z ) Pr ^N 2 (s 1 I z ) Prf 

(hh)pr = £l< V2>*r + A 12; 



(33) 
(34) 
(35) 



the last equation holds for t x s U only. From Eqs. (34) 
and (35), we deduce 



(l-Bi^KA^r-A^, t x >t z . 



(36) 



We symmetrize Eq. (36) by operating on both sides with 
the operator (1 - £ 2 A r 2 ). Because ^ is an instantaneous 
kernel, whereas A 2 is causal, the resulting equation still 
holds for t t ^ t 2 . Also, because A 12 is instantaneous, 
the product i^A r 2 A 12 is zero for t x ^t 2 . We have there- 
fore 



(37) 



for t x a: t 2 . But because Eq. (37) is symmetric, it holds 
irrespective of the time ordering. The effects of pho- 
ton-photon correlation are then given explicitly by Eq. 
(33), (34), and (37) in the form 



(IJz)pr = (1 - £i^i)" l (l - ^* 8 r l Aia, 
( Sl s 2 ) pr = (1 - N&THl - NzB^N^Au. 



(38) 
(39) 



Going back to the unreduced two -point function, we find 
from Eqs. (16) and (39) 

<SiS 2 >, = (1 - - N&Y^az + N X N Z ^ . (40) 



In assessing the performance of the whole system, we 
are finally interested in the statistical fluctuations of the 
path-length error s under the combined influence of at- 
mospheric and photon noise. We use the double bracket 
(( )) to denote an average over both atmospheric and pho- 
ton statistics. We assume that the feedback algorithm 
is chosen to make the mean error 



«s» = (l-A'B)- l <ff> a =0, 



(41) 



which condition will usually be satisfied trivially, be- 
cause the mean atmospheric perturbation (a) a will nor- 
mally be zero. Then the double average applied to Eq. 
(40) gives, by virtue of Eqs. (7), (9), and (29), 

= (1 -N&THl - JWHffi» + A r i Nz^zPA- (42) 

This equation displays the statistical fluctuation of s as 
a sum of two parts, the atmospheric noise term U IZ re- 
duced by negative feedback, and the photon-noise term 
6 12 P Q amplified by the feedback. To recapitulate the 
meaning of the symbols on the right of Eq. (42), A r is the 
feedback kernel, B and P Q are the optical response ker- 
nel and the unperturbed ir radiance that appear in Eq. 
(9), L\ 2 is the atmospheric correlation function defined 
by Eq. (7), and 6 12 is defined by Eq. (30). 

II. OPTIMIZATION OF THE FEEDBACK ALGORITHM 

We take as the criterion for optimizing the system a 
weighted average of a quadratic function of path-length 
errors, 

£=Z Tf&i' **> « s &i> '> s ^> « »> (43) 
where C{k u k 2 ) is an arbitrary symmetric positive-def- 
inite matrix. The value of D is independent of the time 
t at which the averages are taken. We define the system 
to be optimized when D attains its minimum value. In 
the abbreviated notation of Sec. I, we have 

D = ((sCs)), (44) 

where C now denotes the space-time kernel 

C = C{k u h z )${t x -i z ). (45) 

Different choices for C will emphasize different optical 
characteristics of the image. 

As a standard of comparison by which to assess the 
significance of D, we take the quantity 



A,*rf 2 (TrO, 



(46) 



where d is one-tenth of an average optical wavelength. 
This D Q is the value that D would take if each path-length 
error s(k) were an independent random variable with 
variance d z . Because one-tenth of a wavelength is con- 
ventionally considered to be a tolerable error in optical 
systems, the condition 

D< D 0 (47) 

is roughly what is required for the system to give dif- 
fraction-limited images of adequate quality. 

The optical response of the system is sensitive only to 
path-length differences and not to absolute path lengths. 
A path-length error s independent of k produces no 
change of irradiance at the sensors. Formally, by Eq. 
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(9), if u is the constant function 
u{k, t) = l, 



then 



Bu = 0. 



(48) 



(49) 



Because constant path-length errors do not degrade the 
image quality, it is reasonable to make the criterion D 
also insensitive to such errors. We therefore require 
the kernel C to satisfy the condition 

Cw = wC = 0. (50) 

When we substitute the final result Eq. (42) of Sec. I in- 
to Eq. (44), we obtain an expression for D as the trace 
of a product of operators 

D =Tr {C(l - NBY 1 (U+ NI Q N*)(1 - B 'N'T 1 } f (51) 

where the trace indicates a summation over mirror- 
patch indices k and detector indices j and an integration 
over all time coordinates except one. The integration 
over the last time coordinate is omitted because the ex- 
pression in Eq. (43) is independent of time. In Eq. (51), 
B t and N* are the operators adjoint to B and N, defined 
by 

B t {y,x) = B(x,y), N t {x,y)~N{y, x), (52) 
and /„ is the diagonal operator defined by 

with 6 12 given by Eq. (30). Because absolute path 
lengths are immeasurable and only path differences are 
significant, it is often convenient to define the atmospher- 
ic correlation function Z7by 

U(k u k Zy t x - t z ) -.-*<{«(*!, t x )~ aifa, t z )}\ (54) 

rather than by Eq. (7). The expression in Eq. (54), 
unlike Eq. (7), can be directly related to observations 
of atmospheric turbulence. The two definitions of U 
give identical results when substituted into Eq. (51), by 
virtue of Eqs. (49) and (50). It is important not to be 
misled by the minus sign that appears in Eq. (54). In 
spite of the minus sign, U remains a positive-definite 
operator when it operates in the subspace of vectors 
orthogonal to u. 

We now begin the process of minimization of D. Our 
preliminary objective is to prove the following duality 
theorem, which asserts that the feedback kernel N(x, y) 
in an optimized system is proportional to the adjoint of 
the optical-response kernel B(y, x). 

Theorem 1. In an optimized system the feedback ker- 
nel N has the form 

N=KB t /q 1 , (55) 

where K is a causal kernel still to be determined. The 
quantity D defined by Eq„ (43) has the value 

D = Tr {C(l -KT)- 1 {U+KTK')(1 - TK (56) 

with the symmetric kernel T given by 



T = B*Iq 1 B. 



(57) 



posed in Sec. I to ensure that the unperturbed state of 
the system be an equilibrium state. 

Remark 2 e The main problem both in proving and ap- 
plying Theorem 1 is to make sure that the inverse oper- 
ator I 0 _1 in Eqs, (55) and (57) is correctly defined. The 
physical meaning of this Iq 1 is clear; it means that we 
give greater importance in the feedback algorithm to 
photons detected in channels for which the average count- 
ing rate is smaller. The danger in this inverse weight- 
ing of photons is that unphysical singularities may arise 
in channels for which I 0 is zero. We therefore go 
through a careful discussion of definitions to make sure 
that singularities do not arise, either mathematically or 
operationally. 

Proof of Theorem 1. The inverse operator Iq 1 is 
well defined within the subspace L 0 of vectors v(j) for 
which 

v(i) = 0, whenP 0 (i) = 0. (58) 

We define L x to be the space of vectors v of the form 

v^Ba (59) 

with a arbitrary. In words, L t is the linear space of 
possible distributions of ir radiance in the optical de- 
tectors produced in response to all possible atmospher- 
ic perturbations, according to the laws of classical elec- 
tromagnetisnu Statistical fluctuations of the photon 
counts will, in general, allow the real optical response 
to go outside the subspace L x . It is essential to the ar- 
gument that the optical-response kernel B is a classical- 
ly defined quantity that excludes photon-fluctuation ef- 
fects from consideration. 

We prove first that the space L 0 includes L x . The 
classically computed irradiance P(j), which appears in 
the definition (8) of P 0 and B, is a sum of squares of 
wave amplitudes 



po)=X>,o-, S )f, 

p 



(62) 



Remark 1. We shall also prove that the kernel (55) 
automatically satisfies the condition (13) that we im- 



that arise from different parts p of the source when the 
path- length errors are described by the vector s. In Eq. 
(62), we make the assumption that different parts of the 
source radiate incoherently. The amplitudes A p are dif- 
ferentiate functions of the path- length errors s(k\ and 
Eq. (8) gives 

Po(i)-E[^o)] 2 ? ( 63 > 

p 

--zEVi.o)^ V^o" (64) 

From Eqs. (63) and (64), it follows that for each value 
of j, P Q (j) = 0 implies B jk =0 for all k. If, then, v* is 
any vector orthogonal to L 0 , we have 

v' 3 = 0 9 whenP 0 (i)*0; (65) 

therefore 

v'j = 0, when B Jk * 0 for any k, (66) 
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and therefore v' is orthogonal to all vectors in L x ac- 
cording to Eq. (59). This completes the proof that L 0 
includes L 1o 

The inverse Iq 1 being defined in L t on all vectors of 
the form (59), the operator product 

R^tfB (67) 

is defined everywhere and has an adjoint R* . We then 
define the expression on the right-hand side of Eq. (55) 
by the convention 

B'tf-R*, (68) 

which means that the operator (B^q 1 ) is by definition 
zero when it operates in the subspace for which Iq" 1 is 
undefined. The expression (57) for T is also well de- 
fined, and 

R% = B\ R t B= T, (69) 

are ordinary operator equations. 

Let now L z be the linear space of vectors 

v^Ta (70) 

with a arbitrary. Within the space L z , T is a real sym- 
metric matrix with strictly positive eigenvalues; there- 
fore, T has a partial inverse A, such that 

TA = AT=Q (71) 

is the projection operator onto the space L 2 . If v' is any 
vector orthogonal to L z , we have 



(72) 



But Iq 1 is well defined when it operates on the vector 
Bv', and Eq. (72) implies 

(Bv'Y loHBv') - (v'Y TV =0, (73) 

which means that Bv' =0 for all v r orthogonal to L 2 ; that 
is, 

B = BQ-BAT, (74) 

Given any feedback kernel N y we define the related caus- 
al kernel K by 

K = NBA. (75) 

To prove Theorem 1, we have to show that in an opti- 
mized system the relation (75) can be inverted so that N 
is given by Eq. (55). 

If we define N f by 

N=KR* + N\ (76) 

Theorem 1 states that AT' = 0 in an optimized system. 
Equations (74) and (75) give 

NB = NBAT ~KT- KR*B, (77) 
and therefore 

AT'J3 = 0. (78) 

Thus N f is the part of N that operates on vectors orthog- 
onal to Lj. Equation (76) with Eqs. (69) and (78) implies 

NIqN b = KTK f + N'I Q N ft . (79) 

Assembling Eqs, (51), (77), and (79), we find 



D = Tr{C(l - KT)~ l (U+ KTK* + N'I Q N n )(l - TK*)' 1 } 0 (80) 

The iV term in Eq. (76) contributes only to the photon 
noise in Eq. (80) and does not produce any useful feed- 
back. If N f ±0, we can improve the performance of the 
system by setting N' = 0 and dropping the positive-def- 
inite term in N' from Eq. (80). We thus obtain Eq. (55) 
and (56), and so Theorem 1 is proved. 

Proof of Remark 1. We assume that the optical sys- 
tem is geometrically efficient, in the sense that all light 
entering the system is intercepted by optical detectors. 
This does not mean that the detectors are assumed to 
have an optical efficiency of 100%. The optical efficien- 
cy may have any value, provided that it is the same for 
all detectors. The assumption of geometrical efficiency 
means that atmospheric perturbations do not cause any 
light to be lost from the system, although they may 
transfer light flux from one detector to another. For 
any perturbation $(k), the sum of irradiances P(j) given 
by Eq. (8) should be equal to the sum of the unperturbed 
irradiances P Q (j). This means that 

B'w=0, (81) 
where w is the vector of which all components are unity, 

ie(j) = l. (82) 
The definition (53) of / 0 gives 

iherefore Eq. (55) implies 

NP 0 = KB*u: (84) 

The condition (13) is then an immediate consequence of 
Eq. (81), and Remark 1 is proved. 

Our next objective is to determine the time dependence 
of the feedback by optimizing the kernel K. We state the 
result of the optimization again as a formal theorem. 

Theorem 2. In an optimized system with feedback 
kernel N given by Eq. (55), K is the unique causal ker- 
nel that satisfies the nonlinear integral equation 

K+K t -KTK t +U = 0, (85) 

and the minimum value of the quantity D given by Eqs. 
(43) or (56) is 

D^-TriCiK+K*)}. (86) 

Remark 3. It is satisfactory that the optimum kernels 
K and N are independent of the arbitrary matrix C that 
appears in the criterion of optimization D m 

Remark 4. The causality condition to be satisfied by 
Kis 

K = K{k lf kg, t t -t 2 ) = 0 for t x <t z . (87) 
When the trace on the right-hand side is written explicitly, 
Eq. (86) becomes 



(88) 



where #(0+) means the limit of K(t t - t 2 ) as — t z from 
above. Equation (88) shows a direct relation between 
the mean-square path-length errors of an optimized sys- 
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tern and the strength of the instantaneous negative feed- 
back. 

Remark 5. Theorem 2 contains three assertions: 
(a) existence of a causal K that satisfies Eq. (85), (b) 
uniqueness of K, and (c) the statement that the K so de- 
fined minimizes D. We shall prove these assertions in 
turn* 

Remark 6. We already remarked upon the similarity 
between Eq 0 (85) and the Gelfand-Levitan equation. 8 Our 
proof of existence of K is based on the analysis of the in- 
verse scattering problem by Newton and Jost, 9 who gen- 
eralized the Gelfand-Levitan theory to multichannel po- 
tentials. 

Proof of Existence. If W is any kernel, we define P a W 
to be the anticausal part of W. More precisely, 



(89) 



PjMtu h) = W(t u t z ) for t^t 2 , 

P a Wih, <b)=0 for t x >t^ 

Note that P a operates on kernels and not on vectors. Now 
consider the linear integral equation 

G^PJJTG^PJJ ( 9 °) 

as an equation for the unknown anticausal kernel G*. 
This is an equation of Fredholm type (see Riesz and Sz.- 
Nagy 10 ). The operator P a UT is completely continuous 
when applied to G*. Therefore the Fredholm alternative 
holds: either the inhomogeneous Eq. (90) has a solution, 
or the corresponding homogeneous equation 



F+P n UTF = 0 



(91) . 



has a nonzero solution F. But Eq. (91) implies for an 
anticausal F 

Tr[F*TF] + Tr [F'TUTF] = 0 (92) 

and, since both T and U are positive definite, Eq. (92) 
implies TF = 0. Then Eq. (91) gives F = Q, so we have 
proved the existence of G* satisfying Eq. (90). We de- 
fine the kernel K by 

K^Gt+UTGt-U, (93) 
so that Eq. (90) becomes 

P a K=0. (»4) 

Thus K is a causal kernel; we now have only to prove 
that it satisfies Eq. (85). 



The equation adjoint to (93) is 
K*=G+GTU-U. 



(95) 



We multiply Eq. (93) on the left-hand side by (1 - GT), 
multiply Eq. (95) on the right-hand side by (1 - TG% 
and subtract. The result is 



G + K- GTK= G f +K* —K t TG t m 



(96) 



But the left-hand side of Eq. (96) is causal, whereas the 
right-hand side is anticausal. Therefore both sides are 
zero, and we have 



(1 - GT)(1 - KT) = (1 - TG)(1 - TK) = 1, 
G= - (1 - KT)'% GT)~ X G, 



(97) 
(98) 



the inverses here being well defined by virtue of Eq. 
(97). When we substitute G from Eq. (98) into (95), the 
result is Eq. (85). This completes the proof of the ex- 
istence part of Theorem 2. 

Proof of Uniqueness. Since T and U are positive def- 
inite, the operator (1 + UT) has an inverse. Now if K is 
any solution of Eq. (85), we have 



{l-KT)(l-K t T) = l + UT , 



(99) 



and so the inverse (1 — KT)" 1 also exists. If K x and K z 
are any two solutions of Eq. (85), we find by equating 
the two expressions for U, 

(1 - K^YHK, - K 2 ) = (Kt - K[)(l - TKi)' 1 . 



(100) 



If K x and K z are both causal, the left-hand side of Eq. 
(100) is causal and the right-hand side is anticausal, so 
that both sides must be zero. This proves the unique- 
ness of K. 

Proof that D is a minimum. Let K be the solution of 
Eq. (85) and let K* be any other causal kernel for which 
the value D r of D is defined by Eq. (56). We define a 
causal kernel H by 

H = {l-K'Ty l (K'-K), (101) 
so that 

(1 - K'TY 1 = (1 +HT)(1 - GT) , (102) 

(1-iC'T)" 1 K' = H-(1+HT)G . (103) 

Using Eq. (102) and (103), we eliminate K' from Eq. (56), 
and express D' in terms of G and H. We eliminate U 
from D f by using the equation 

(104) 



(1 - GT)U{1 - TG*) = G + G f - GTG< , 



which is obtained by substituting for K from Eq. (98) into 
Eq. (93). The resulting formula for If is 

D r = Tr{C(G +G f +HTG + G t TH t +HTH*)} . (105) 

The purpose of the substitution Eq. (101) was to reduce 
D' to a quadratic form in the unknown H. Now 

Tr {CHTG}= Tr {CG'Ttf*} = 0 , (106) 

because both G and H are causal kernels while C and T 
are instantaneous. Therefore D' depends on H only 
through the quadratic term which is positive definite. So 
we have proved that D' attains its minimum value when 
H= 0 and K' =K. The value of D' at the minimum is 

2) = Tr{C(G + G'} . (107) 

We have also 

Tr{CGTtf}=Tr{Ciir'rG'}==0 , (108) 

because G and K are both causal. Because both sides of 
Eq. (96) are zero, the value of D given by Eq. (107) co- 
incides with Eq. (86). This completes the proof of 
Theorem 2. 

III. APPLICATION TO PARTICULAR SYSTEMS 

The theory in Sees. I and n can be applied to active 
optical systems of many kinds. We have carried the 
analysis through in detail for three special cases, (1) a 
simple stellar interferometer of the type developed by 
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Currie, u (2) a conventional telescope with detectors in 
the image plane, as described by Muller and Buffington 2 
and by Bridges et al., z and (3) a system that uses an 
array of interferometers, as described by Hardy et at. 4 
and by Dicke. 12 We take the criterion of successful op- 
eration of the system to be Eq. (47), with D and D 0 de- 
fined by Eqs. (43) and (46). According to this criterion, 
a successful system should give angular resolution equal 
to the diffraction limit. In each case, we suppose the 
feedback algorithm optimized according to the prescrip- 
tions of Sec. n. The following is a brief summary of the 
results. 

(1) For the simple interferometer, we have as the con- 
dition for successful operation 

(Pp z rfre>a z r (109) 

where <j> is the incoming photon flux (photon number per 
unit area per unit time), p is the diameter of each of the 
two receiver apertures (assumed small enough so that 
atmospheric phase differences within each patch are un- 
important), rj is the fraction of the light contributing to 
optical structure at the desired resolution (equal to the 
fractional modulation of the observed interference pat- 
tern), r is the average autocorrelation time and a the 
root-mean-square amplitude, in radians, of the atmo- 
spheric phase difference between the two interfering 
beams, and e is the fractional efficiency of the detectors. 

(2) For the system with image -plane detection, the 
. criterion for success is 

4>p z T?Te>{A/p z ) , (110) 

with 0, n, c defined as before, whereas p is now the 
seeing-patch diameter (diameter of the patches on the 
mirror that are independently controlled), r is the auto- 
correlation time of the phase differences between neigh- 
boring patches, and A is the total area of the mirror. 

(3) For a system with detection by multiple interfer- 
ometers that monitor phase differences over all of 
the receiving area, the criterion for success is 

<t>p z r} z re>l i (111) 

with the symbols defined as in case (2). In each of Eqs. 
(109), (110), and (111), a numerical coefficient of the 
order of unity has been dropped, the precise value of 
the coefficient depending on details of the system design 
and of the atmospheric fluctuation spectrum. The quan- 
tity that appears on the left-hand side in each case is the 
effective number of photons that are detected per reso- 
lution element of the system in space and time. 

The criterion (110) for image -plane detection has the 
aperture area on the right-hand side. This means that 
image-plane detection is undesirable for large systems. 
The inefficiency of image-plane detection arises from 
the mixing of photons from all parts of the mirror at 
each of the detectors. When photons are mixed, infor- 
mation is lost that cannot be recovered in subsequent 



data processing. The condition (111) shows that with 
interferometric detection the limiting photon flux is in- 
dependent of the size of the system. If an object is 
bright enough to fulfill condition (111), it can be observed 
with angular resolution limited only by the size of the 
instruments that it is feasible to build. It is likely that, 
in optical astronomy as in radio astronomy, the very- 
high-resolution instruments of the future will be inter- 
ferometric arrays rather than single dishes. 

If we substitute into Eq. (Ill) numerical values appro- 
priate to conditions of good seeing at ground -based ob- 
servatories, namely 

£ = 10 cm, t-10" 2 s , (112) 

with a source possessing sharply defined optical struc- 
ture (77 = 1), and assume a detection efficiency e = 0. 1, 
we find the limiting source brightness to be 

0 = 10 ernes' 1 . (113) 

This means that objects brighter than magnitude 14 
should be observable with arbitrarily high resolution, 
under conditions of good seeing. The numerical esti- 
mate is of course very rough, and is intended only as a 
guide to the planning of experimental systems, rather 
than as a prediction of their performance. 
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Least-square fitting a wave-front distortion estimate to an array of 

phase-difference measurements 

David L. Fried 
Optical Science Consultants, P. O. Box 446, Placentia, California 92670 
(Received 17 July 1976) 

The problem of fitting a wave-front distortion estimate to a (single -instant) set of phase-difference 
measurements has been formulated as an unweighted least-square problem. The least-square equations have 
been developed as a set of simultaneous equations for a square array of phase-difference sensors, with phase 
estimates at the corner of each measurement element. (This corresponds to the standard Hartmann 
configuration and to one version of a shearing interferometer of a predetection compensation wave-front 
sensor.) The noise dependence in the solution of the simultaneous equations is found to be expressible in terms 
of the solution to a particular version of the measurement inputs to the simultaneous equation, a sort of 
"Green's-function" solution. The noise version of the simultaneous equations is solved using relaxation 
techniques for array sizes from 4 X 4 to 40 X 40 phase estimation points, and the mean-square wave-front error 
calculated as a function of the mean-square phase-difference measurement error. It is found that the results 
can be approximated within a fraction of a percent accuracy by <(8<t>) 2 > = 0.6558[1 + 0.2444 InCN 1 )]^, 
where <(54>) 2 > is the mean-square error (rad 2 ) in the estimation of the wave-front distortion, for a square 
array consisting of N 2 square subaperture elements over which two phase-difference measurements are 
made — one phase difference across the x dimension and the other difference across the y dimension. Here 
cr^ is the mean-square error (rad 2 ) in each phase-difference measurement. 



INTRODUCTION 

In this paper we shall address the question of the per- 
formance of a processor which, given a set of phase- 
difference measurements across a large aperture, has 
to formulate a corresponding estimate of the distorted 
wave front over that aperture. This processor is the 
key to certain types of adaptive optics systems, par- 
ticularly those that utilize a pseudo- Hartmann or a 
shearing- interferometer- type sensor. 1-8 

We measure the performance of the processor in 
terms of the rms error in the estimated wave front, 
((a*) 2 ). We assume that the mean-square error in the 
input phase-difference measurements, a^, is known 
and will see that it controls the error in the processor 
output. We further assume that since there are more 
phase-difference measurement inputs than points at 
which the processor is to estimate the phase difference 
(nearly twice as many), a least-square error estima- 
tion procedure will be used. 

DEFINITION OF ELEMENTS ON THE APERTURE 
PLANE 

For the purposes of our analysis, we shall consider 
an aperture in the form of a square array of square 
regions, arranged like a checkerboard. We shall re- 
strict our attention to arrays consisting of an odd num- 
ber of elements, i.e., the array size will be 

tf 2 =(2n+t) 2 . (1) 

We consider the element to be a region across which a 
phase difference is estimated. For the pseudo-Hart- 
mann measurement technique, the measured phase 
difference between opposite sides of the element would 
correspond to the measured angle-of-arrival times the 
element dimension, converted from length into radians 
of phase. We would write 

<P x = kd9 x , (2a) 

<P^kdO y , (2b) 



where k = 2ttA is the optical wave number, d is the 
length of a side of the element, 9 X and 6 y are the mea- 
sured angle-of-arrival components, and 0 r and <p v are 
the phase differences between the two sides of the 
square element. If a shearing interferometer is used 
to measure the phase difference, then it is convenient 
to consider the shear distance to be 5. Then 

* r =A<Mrf/fi), (3) 

(4) 

where A0 r and A# y are the measured phase differences 
for shear of distance 6 in the x and y directions, re- 
spectively. 

It is not our intention to take up here the question of 
the relative merits of each of these phase-difference 
measurement options. There are considerations con- 
cerned with photon shot noise effects and the fact that 
the measurements are only approximately the phase 
difference across the aperture element rather than the 
true phase difference between the two sides as a result 
of spatial quantization in our measurements. At this 
point we simply wish to know the rms error in the esti- 
mated instantaneous wave-front distortion, assuming 
that we have meaningful measurements of the phase dif- 
ference, <p x and <p y9 between the two sides of the aper- 
ture element. 

We shall define the instantaneous wave- front distor- 
tion by phase estimates at each of the four corners of 
each aperture element. It is convenient to denote the 
element by the notation (i 9 j) with -n^i**n and - « < j 
< n. We shall let <p Xtttj and <p^ Uj denote the phase dif- 
ference between the two sides of the (:,;")th aperture 
element. We shall use the notation to denote the 
estimated phase (matching the instantaneous wave- 
front distortion) at the corner common to the (*=£+£, 
j = (*=/>+!, ; = {i = p-i, j = q+i), and 

(i = p - 2 , j = q - i) aperture elements. We let p and q 
run through the values p = -n- i, -n+i, -n+|, ... > 

i_ _1 3 3 1 1 . . -ii 

— 2, 2 , 2, ... , »-?, n — 2 , «+2, ana similarly, 
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0. J>th aperture 
element 



FIG. 1- Representatioa of the orientation and location of phase 
differences and phase estimates for the (i,;)th aperture ele- 
ment. 



«- 5, «+2. The relationship between 0 rW< y , 



and the various for the (z,;)th aperture element is 
indicated in Fig. 1. 

With the quantities of concern to us in our problem 
thus defined, we now turn our attention to the formula- 
tion of a method of calculating the § Ptq terms that make 
up the instantaneous wave-front distortion estimate, 
given the instantaneous measurements of the phase dif- 
ferences, <p XtitJ and <P ywltJ . This is treated in the next 
section. After that we shall turn our attention to mat- 
ters related to the estimation of the rms error to be 
expected in * hq9 given the rms error in 0 Xt(t/ and in 



LEAST-SQUARE ESTIMATION OF <P 

We note that there are two phase-difference mea- 
surements, * Xflti and <f>y, UJ , for the (ij)th aperture 
element, so that for ;V 2 aperture elements there will 
be 2iV 2 measurements in each instant. We also note 
that for ,V 2 aperture elements there are (N+ l) 2 values 
of ^ Ptq to be estimated. The $ Ptq are degrees-of- free- 
dom to be adjusted to match the constraints provided 
by the <p XfUf and <£ yWfi measurements. We see that for 
-V^3 there are more constraints than adjustable de- 
grees-of- freedom, and so a least- square fit to the data 
is called for. 

The phase differences across the (ij)th aperture ele- 
ment are related to the phase estimates by the equa- 
tions 



*x, I W = 2 (* f ♦d/2), }+a/Z) + $ { +U/2), /- (1/2)) 
^Vff./^ 2(*<*U/2)./*U/2)+$f-<l/2)./*(l/2>) 



(5a) 



J 



- 2(*f*(l/2) t /-(i/2)+^f-(l/2) t i-(l/2)) , (5b) 

where 4> XtUj and $ y%Uj are defined by Eqs. (5a) and (5b), 
with the $ Pfq chosen so jis to minimize the discrepancy 
between the i Xfit3 and $ ytltJ as compared to the <P Xritj 
and $y,i,j- Tnis discrepancy can be defined in un- 
weighted least- square error terms as 



it 1 



(6) 



We wish to choose the $ Ptq so as to minimize the 
squared error, A. 

It should be noted that if we were only interested in 
short-exposure imagery, for which the average tilt 
across the aperture has no effect, then Eq. (6) would 
have to be modified to reflect this. This would be ac- 
complished by subtracting the average over of 
0r,*\y* -0 r , fM . from each O^,,,,- <P s , itJ ) term before 
squaring, and similarly subtracting from the 
- 0 yt * w ) terms before squaring. Since, however, we 
are interested in long-exposure resolution, failure to 
consider an average tilt error would be erroneous. We 
therefore have left the average tilt error in our expres- 
sion for A as formulated in Eq. (6). 

Our selection of the estimated phase values is based 
on the requirement that A be minimized, Thus we 
write 



3A 

9<K, 



(7) 



for all values of (p, q), thus establishing a set of (AT+ I) 2 
simultaneous equations to be solved for the (.V+ I) 2 un- 
known values, If we substitute Eq. (6) into Eq. 
(7), we get for all (p,q) 

(8) 

Making use of Eqs. (5a) and (5b), we see that 



3<f> 



- Hp, i - $)5(q,j + i) -5(p,i- iWq,j - k) , 



(9a) 



-Hp, i+i)5(q,j- i) -6(p,i -i)6(q,j-i). 



The delta functions here are intended to denote 
Kronecker deltas, according to which 



10 if a*b. 



(9b) 



(10) 



If we substitute Eqs. (9a) and (9b) into Eq. (8), and ex- 
ploit the property of the Kronecker deltas to allow us 
to carry out the (ij) summation, we obtain for all (p, q) 



-U/2W(l/2) - 0x,,-<1/2>.,-(1/2)) + W>r.,-(l/2WU/2) - <? y , ,-a/2WU/2>) " W> r . *U/2WU/2) - *x.Ml/2 WU/2)) 
~ (0 *-*<"2W«/2>-<^^ 

" ( *y.#-<l/2).,*(l/2) - *»#-U/»,«»U/2)) - (<f\.*U/2).,*U/2) - *y,*(l/2 WU/2)) = 0 . 



(ID 
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If we substitute Eqs. (5a) and (5b) into Eq. (11) and coLlect the various $ terms, we get for all (p, q) for which 
p*±(n + i) and q*±(n+i) 

*^-l f9 -r *A-W«*l - ^0*l.<rl " $ M,«*L = ^Ml/2).ra/2)+^^-«l/2),«»a/2) ~ 0r.JH(l/2),«-(l/2) ~ 0r, U/2 WU/2) 

+ ^y,J>-a/2),«r-(l/2)+ ^y f ^d/2),<r-(l/2)- 0y, *-U/2), «+U/2) ~ 0y, >+U/2) f «*U/2) - (12) 



For the various cases in which p and/or q equal ± («+ 2), 
certain terms in Eq. (12) would be missing since they 
pertain to aperture elements that are not part of the 
(2n + l) 2 elements making up the total array. There 
are eight variant forms of Eq. (12) governing the cases 
for which p and/or q equal ± (n+ £). These are most 
compactly represented if we define the functions f(i,j) 
and g(i,j) as 



1, 


if |f| *zn+i and \j\ 






otherwise ; 






if |t| <n+i and \j\ 




1, 


if |i|=n+| and |j| 




2, 


otherwise . 





(13) 



(14) 



Making use of these functions, we can write the general version of Eq. (12), which is valid for all values of {p, q) as 

I, <7- i)(0x.^U/2K«-a/2)+ ( f > 3P.^-(l/2),«-(l/2)) +/(/>- i ?+^)(^r,j>-(l/2)t<i*(l/2)- 0»»-<l/2).«+U/2)) 
+/(/> + St <7~ *r.*U/2),«-U/2> + 0y.^(l/3).^a/2)) +/(/>+ 2, * x , ^(1/2), ^U/2) - #y,*U/2>,«*U/2>) . (15) 

Equation (15) represents (2n+2) 2 simultaneous equations giving the values of the (2n + 2) 2 values of $0,,, based on 
the 2(2rc + I) 2 values of the measured phase differences <t> XtUj and <b ytUJ . These simultaneous equations represent 
a least-square fit of the instantaneous estimated wave-front distortion to the phase-difference measurements. 

With the least- square solution formulated in terms of Eq D (15), we are now ready to turn our attention to the 
evaluation of the mean- square instantaneous wave-front distortion error to be expected as a consequence of some 
rms phase-difference measurement error". We take up the formulation of this problem in the next section. 



rms WAVE-FRONT ERROR— FORMULATION 

We start our evaluation by noting that Eq. (15) repre- 
sents a set of {2n + 2) 2 linear simultaneous equations. 
The fact that the equations are linear means that if 

are a set of (2n + 2) z solutions associated with 
the set {(0 Xr f,Pii W>y,i,/)i} of 2(2«+ I) 2 measurements, 
and if {(^^J is the set of solutions corresponding to 
the measurement set{(0 r?{ } ) z ; (<t> VfitJ )^, then the set 
of solutions {(3\ c )i+ (fc*,)^ corresponds to the mea- 
surement set {(*». ff ,)!+(*„ i, (** , fi )i+ tt>»i t/ )J. 
Recognizing the implication of this linearity, we per- 
form the following decomposition of the measurements 
and calculated instantaneous wave-front distortion: 



(16a) 
(16b) 



where <t>^ itj and <t>^ Uj are the correct values that the 
measurements should yield (if there were no noise), 
and 6# Xtiti and 54> y>iti are the noise contributions to the 
measurements. We write the calculated instantaneous 
wave-front distortion 



(17) 



where 4>^ a are the correct values of the instantaneous 



wave-front distortion, which would be calculated by 
using the <P TtitJ and 4>^ {tjf so that 5^ A?C is the error in 
the calculated instantaneous wave-front distortion due 
to the measurement noise . It follows from the linear- 
ity of the simultaneous equations that {6$ N J is the 
solution of Eq. (15) if {&<P x , l9j ; are taken as the 

measurements to be input into Eq. (15). 



! 

The implication of the linearity of the simultaneous 
equations can be pushed even further. If we take all 
{($p, q ) x ,i,j} to be the solution to Eq. (15) when the 
2(2n+ l) 2 measurements are replaced by zero, except 
for the <b x . ltJ measurement which is replaced by unity, 
and if {($ Nc ) y ,i ? j} is the solution when all the measure- 
ments are zero except the <t><,, us measurement, then it 
follows from the linearity of the equations that 

5**,= E 6 ^ i .i(*N,)x.i.y+ 6 ^yW.i^^,)yW.i. d 8 ) 



The mean- square wave-front distortion can be writ- 
ten 



««*)*>= (2n+ l)- 2 <?K(6**,) 2 >, 



(19) 



where we have weighted the (p, q) locations on the in- 
terior of the array covering the aperture at unity, 
those at the four corners at one-quarter, and those on 
the edge but not at the corner at one-half. With this 
weighting to appropriately discount (p, q) boundary 
points, the factor (2rc+ l)" 2 provides the correct nor- 
malization for the (2rc+2) 2 term summation. 

The evaluation of ((6$ N ,) 2 ) follows directly from Eq. 
(18). We shall make use of the statistical independence 
of the measurement noise values, and the fact that each 
measurement has the same rms phase-difference error, 



. We write this in the form 
( 5 *y,^i 6 *y.i'.r> = cr « 5 ^ £ ' )6 0',J / ), 



(20a) 
(20b) 
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,6^,j^r,,)-0, CMC) 
re 5te & ) is the Kroneck:er del ta, defined in Eq, 
*0) ^ w ' e subsfcitute Eq - (18) in the evaluation of 

) 2 )> we can cast the result in the form of a 
adruple sum on i, but by making use of Eqs. 

q ?0a) (20b), and (20c) this can be reduced to a double 
s'umon**,;. Thus we obtain 

((5^ N /> = C7^g{[(^,) T . t -./J 2 ^(*NAi./}. (2D 

Vow if we substitute Eq. (21) into Eq. (19) and inter- 
change the order of the (p, q) summation and the 
summation, we get for the mean- square instantaneous 
w ave-front distortion error 

<(*♦)*>= (2n + ir 2 Z (°* E ±*P> 9){[{^Xu S f 

+ [(*j»«.,n)- (22) 

At this point as a matter of practical simplification, we 
introduce the assumption that the sum over (/?, q) has a 
value which is very nearly independent of (i,j) 9 i, e. , 
we assume that the wave-front distortion mean- square 
value over the aperture due to a single measurement 
error is independent of where the error is located, and 
also of whether it is a <P XtitJ error or a error. 
This allows us to replace [{$p,X.uj\ Z and [(**«)**,/]* 
in Eq, (22) with the corresponding values for i = 0, j-0, 
with no need to distinguish between subscript x and 
subscript y. We write 

»*Joo« (**•>*!., fori = 0, ; = 0, (23) 
and now can cast Eq. (22) in the form 

((o*) 2 >= (2n + I)" 2 Z *>[(**,)»]*) • (24) 

We note that the quantity in the large parentheses is 
independent of i 9 j 9 so that when the sum is carried 
out, it simply gives a factor of (2n + l) 2 (since i and; 
each run over the 2n+ 1 values from -n to Thus 
Eq. (24) is reduced to the form 



where 



(25) 



(26) 



The problem of evaluation of the mean- square instan- 
taneous wave-front distortion error ((5$) 2 ) associated 
with an rms phase difference measurement error a M 
is seen to reduce to the evaluation of 31. Results of this 
evaluation for various size aperture arrays are pre- 
sented in the next section. 

rms WAVE-FRONT ERROR—EVALUATION 

The evaluation of the rms instantaneous wave-front 
distortion error due to phase-difference measurement 
errors requires the evaluation of 31, and as we see 
from Eq. (26), this is basically the problem of evaluat- 
ing the set {($ Na ) 0 o}. The evaluation of the (2n + 2) z 
terms in the set {(* Na ) 0Q } requires the solution of Eq. 



(15) with all the <p XtitJ and <P v , Uj replaced by zero, ex- 
cept <£ s ,q,o> which we replace by unity. For the four 
cases defined by p = ± i , q = ±j, the right-hand side of 
Eq. (15) does not vanish. Otherwise it does. Thus Eq. 
(15) reduces to the five specific cases of 

(*i/2,i/2)oo = [5 r (i, a)]" l {/(- a, i)(*-i/ 2 ,.i/ 2 )oo 

+/(- 2, i)($-l/2, 3/2)a0+/(t, - 3/2,-1/2)00 

+/(i?)(*3/2,3/ 2 )oof l}, (27a) 

(* 1/2.-1/2)00 = [g(k, - i)]~ l {/(- i - !)(*-i/ 2 ,. 3 / 2 )oo 

+/(- 2, 2, 1/2)00 +/(t, - 3/2. -3/2)00 

+ (i*)(*a/i t i/Joo+l}, (27b) 

(*-l/l B l/t)00- i *>l" l {/(- i - 2-)(*-3/2.-l/ 2 )00 

+/(- 1, i)(^-3/2.3/ 2 )00+/(i - 5)(*l/2.-l/z)oQ 

+/(il)(* 1 / 2 .3/ 2 )oo-l}, (27c) 

(*-l/2.-l/ 2 )00=[^(- i-i)]-H/(-i-i)(^3/ 2 .-3/ 2 )00 

+/(" f, 4)(*-3/2. 1/2)00 - a)(* 1/2.-3/2)00 

+/(i*)(*i/a.i/i)oo-l}. (27d) 

and for (p, <?) * (± i, ± i), 

(**«)oo = [giP, <?)Y l {f(P~ 1, tf- D(*,-i. 9 -i)oo 

1, ? + l)(*,. Ut Joo+/(*+ 1, »(* *i it .i)oo 
+/(P+U+l)(WiU. (27e) 
Equations (27a)- (27e) are in a form that appears rea- 
sonably well suited to solution by relaxation techniques. 
We have found, however, that application of straight- 
forward relaxation methods does not seem to yield a 
satisfactory result. We suspect that this is related to 
one or both of the following facts about this set of equa- 
tions. (1) The coefficients of these simultaneous equa- 
tions are such that if {($ Nff )oa} is a solution set, then so 
is {(Sfcjoo+C} where C is any constant. (2) There are 
actually two sets of equations that can be solved sepa- 
rately. The two sets consist of those ^ Ptq for which 
p+ q is odd and those for which it is even. Each set 
can be solved independently, and then different con- 
stants, and C„ en , can be added to the sets of solu- 
tion values. What this means is that the equations as 
written do not have any continuity between horizontally 
or vertically displaced adjacent elements, only between 
diagonally displaced adjacent elements. 

To avoid these problems, to get the equations in a 
form where relaxation methods could be successfully 
applied, and incidentally to significantly reduce the size 
of the computational task, we have taken note of the ex- 
pected symmetry of the solution to set up a reduced- 
size problem, 

Our approach is to argue that by virtue of the sym- 
metry of the problem, 

• (28) 

Also, if we restrict attention to a solution in which the 
average value of <£ Nfl is zero (this restricts our choice 
of C, the additive constant), then 

(29) 
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TABLE I. Mean-square noise coefficient for fitting a wave- 
front distortion estimate to a set of noisy phase-difference 
measurements. The significance of the noise coefficient, ai t 
is expressed by Eq. (25). 



K 


31 


iV 2a 


(2/0 2b 


2 


0.4900 


9 


16 


3 


0.5832 


25 


36 


4 


0.6398 


49 


64 


5 


0.6810 


81 


100 


6 


0.7136 


121 


144 


7 


0.7405 


169 


196 


8 


0. 7634 


225 


256 


9 


0.7834 


289 


324 


10 


0. 8012 


361 


400 


15 


0.8684 


841 


900 


20 


0.9151 


1521 


1600 



z N z = (2K-l) 2 is the number of measurement elements in the 
array. There are 2N l phase-difference measurements made. 

\2K) 1 is the number of points at which the wave-front distor- 
tion is calculated. 



This means that where previously we had to consider 
(2«+ 2) 2 unknowns for a (2«+ l) 2 set of phase- difference 
measurement elements, we now only have to consider 
(« + l) 2 unknowns. This greatly reduces the size of our 
computational load, We see now that we will have 
basically different equations for the case of p = j and/or 
q = j than for other cases, since a quantity such as 
^V-wi 011 fche right-hand side of Eq. (21 e) would be re- 
placed by- for/? =4, or $ pmUq ior <?=i, with corre- 
sponding corrections for Eqs. (27a)-(27d). This, inci- 
dentally, would eliminate the possibility of separating 
the problem into two independent sets of simultaneous 
equations, thus insuring continuity in our solution set. 

Making use of Eqs. (28) and (29) to modify Eq. (27a), 
we get ' 

(*i/2,i/2)oo=[£"(i S>]"H-/(i s)(*i/2,i/ 2 )oo 

-/(I, l)(*l/2 t 3/2)00+/(i £)($3/2,l/ 2 )00 

i)(* 3 /2.3/ 2 )oa+l}. (30) 
Equation (30) can be further simplified to take the form 

(W/2>oo=[<r(i,i) i)(*i /2 , 3/2 )oQ 

iH$3/2,l/2)00+/(i I )(*3/2,3/2)00+ U • 

(31a) 

We note that there is no need to seek modified versions 
of Eqs. (27b), (27c), or (27d) since the same informa- 
tion can be obtained from the solution of E-\. (31a), by 
use of Eqs. (28) and (29). Equation (27e) is subdivided 
into the following three forms accordingly as p = i or 
q=2 f or neither equals one-half. We get 

(*i/&a)oo=[*(i ?)]"H-/(i q- D(*i,z ie .i)oo 

-/(i q+ D(* 1/2 ,^ x ) M +/(i q- 1)^3/2,^1)00 
+/(i^+l)(*3/2,«*i)oo} (<7>i), (31b) 

(*M/8>oo= [*(/>, i)Y l {f(p- 1, i)(*p- U v 2 )oa 

+AP-1, t)(*,-i.3/ 2 W/(/>+ 1, i)(*^ U i/«)oo 
+AP+ l,D(*^i t 3/ 2 )oo} (/>>i), (31c) 



and,, of course, 

(*p.Jw=l*(P, q)]~ l {Ap-l,q-l)(* P -u<-i)oo 

+Ap - 1, <7 + i)(*^-Utrt)oo+/(* + If <?~ U(**l.«-i)oo 

1, <? + D(3w +l )oo} (/> > I and <7 > 1) . 

(3 Id) 

In preparing a computer program to solve for ($^) Q0 
and evaluate 31, we have considered a system with 
(2K- l) 2 measurement elements and (2K) Z positions at 
which the phase estimates had to be formed, We there- 
fore had to solve for K z unknowns. [Note: Hereafter 
we use K to denote (w+1), i.e., K=n+l.] The method 
of solution is a straightforward relaxation procedure 
with a full new solution set obtained in each iteration 
from the previous solution set before any of the pre- 
vious solution set elements are replaced. We have car- 
ried out the solution process for K = 2, 3, 4, , . . , 10 and 
for K= 15 and 20. To be sure that our solution had 
converged, the program went through olK iterations 
where a = 10 for K= 2, 3, 4, 5, and 6, a = 20 for K= 7, 8, 
9, and 10, and a = 30 for if = 15 and 20. The program 
calculated and printed out the current estimate of 01 
every tenth iteration. We were able to see that the 
value of 0t had stabilized to better than four significant 
figures in all cases well before we reached the aKth 
iteration, After the aKth iteration, the program would 
print out the solution set, or rather the positive p, 
positive q quadrant of the set, Then it would go through 
one more iteration, and print out the next value of 01 
and the updated solution set. In no case did we find 
any suggestion in examining any of the printouts that the 
iteration procedure had not fully converged, at least to 
within the limits of accuracy imposed by the six-sig- 
nificant-figure nature of the computations, The results 
obtained appear to be quite reliable for our purposes. 

DISCUSSION OF RESULTS 

The values of 01 calculated by the program are listed 
in Table I, If we examine the data in Table I, we find 
that it can be quite accurately estimated by the function 

(32) 



01 = 0. 1603 ln(tf- k) + 0, 4390 (K > 2) 



We recall that we are dealing with a system in which 
there are N z measurement elements (and 2N Z phase- 
difference measurements possible) and (N+ I) 2 points 
at which the phase is to be estimated, Obviously, then, 
we can equate if with the measurement array size by 
the formula 



N 2 =(2K-l) 2 . 



(33) 



If we make use of Eq. (33) to eliminate K in Eq. (32), 
we get 

01 * 0. 08015 IntV 2 ) + 0, 3279 . (34) 

If we substitute this result into Eq. (25), we see that 
the mean-square phase error in estimating the wave- 
front distortion over a square aperture of AT 2 measure- 
ment cells will be 

((S*) 2 ) = 0. 6558[1 + 0. 2444 ln(tf 2 )]o^ , (35) 

where a 2 ^ is the mean-square error expected in mea- 
suring the phase difference in one direction across a 
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te measurement ceil. (It is interesting to note that 
S H\ery simple equation <(5$) 2 > = f[l + * In(,V 2 )]cr^ 

Ids results which appear to be in error by no more 
y ^ about 3%. ) As can be seen from Eq. (35), the 
1 an-square error in the phase estimate is virtually 
^dependent of the size of the array. Apparently the 
0 ise effects are very much a local matter and the 

ence or absence of boundaries to the array in the 
P icinity of each noisy measurement has very little ef- 
fect on the way it introduces error into the wave-front 
estimate. 

ijt should be noted that while there has been an extensive body 
0 £ work, and some previously published resuLts pertaining to 
adaptive optics design and wave-front distortion sensing, a 
few examples of which we list in the references below, al- 
most none of this work considered the approach of measur- 
ing an array of phase differences and thus had to consider 
the data fitting problem that we consider here. 

X. R. O'Meara, U. S. Patent 3,764,213, "Return Wave 
Phase Control Adaptive Array," (9 October 1973). 



3 J. W. Hardy, J. Feinleib, and J. Wyant, ,c Real Time Phase 
Correction of Optical Imaging Systems," in the Digest of 
Technical Papers of the OSA Topical Meeting on Optical Propa- 
gation Through Turbulence, 9-11 July 1974, Univ. of Colo- 
rado, Boulder, Colo. (This paper considers use of a shear- 
ing interferometer to measure phase differences, but there 
is no discussion of the wave-front fitting problem. ) 

4 W. T. Catney, C. L. Hayes, W. C. Davis, and V. F. Piz- 
zurrd, "Compensation for Atmospheric Phase Effects At 
10.6 fi," Appl. Opt. 9, 701 (1970). 

5 R. A. Muller and A. Buffington, "Real-time correction of 
atmospherically degraded telescope images through image 
sharpening," J. Opt. Soc. Am. 64, 1200-1210 (1974). 

6 C. A. Primmerman and D. G. Fouche, "Thermal-Blooming 
Compensation: Experimental Observation Using a Deforma- 
ble Mirror System," AppL Opt. 15. 900 (1976). 

7 W. B. Bridges, P. T. Brunner, 5. P. Lazzara, T. A. 
Nussmeier, T. R. O'Meara, J. A. Sanguinet, and W. P. 
Brown, Jr., "Coherent Optical Adaptive Techniques," Appl. 
Opt. 13, 291 (1974). 

3 J. E. Pearson, "Atmospheric Turbulence Compensation 
Using Coherent Optical Adaptive Techniques, " Appl, Opt. 
15, 662 (1976). 
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APPENDIX G 



Reprinted from Journal of the Optical Society of America, Vol. 67(3), pp. 375-378 

(March 1977). 



Wave-front reconstruction for compensated imaging 

Richard H. Hudgin 
Itek Corporation, Lexington, Massachusetts 02137 
(Received 13 August 1976) 

A critical component in a compensated imaging (CI) system is the wave-front sensor which measures the 
residual distortion of the wave front after reflecting off the active mirror. The sensor produces estimates of 
wave-front slopes or phase difference across the aperture. For many applications, the phase differences or 
slopes are not the most convenient form of data for processing or control, and they must be converted to 
absolute wave-front phases. This paper analyzes the conversion from phase differences to phases and derives 
the optimal linear estimator in terms of least noise propagation. Some remarks concerning hardware 
implementation are also made. 



I. INTRODUCTION 

A critical component in a compensated imaging (CI) 
system is the wave-front sensor which measures the 
residual distortion of the wave front after reflecting off 
the active mirror. The sensor produces estimates of 
wave-front slopes or phase difference across the aper- 
ture. For many applications, the phase differences or 
slopes are not the most convenient form of data for 
processing or control, and they must be converted to 
absolute wave-front phases 0 This paper analyzes the 
conversion from phase differences to phases and de- 
rives the optimal linear estimator in terms of least 
noise propagation. Some remarks concerning hard- 
ware implementation are also made. 1-10 

II. SENSOR MODEL 

For a wave-front sensor divided into many subaper- 
tures, a set of noisy phase-difference estimates are 
produced which measure the phase differences between 
adjacent subapertures. For a two-dimensional array of 
subapertures, there are four phase differences connect- 
ing one subaperture to its neighbors. Let <p jk be the 
two-dimensional array of subaperture phases where j, 



k = 1, . . . , N 9 and let the connecting phase difference be 
defined by 11 

s )k =< Pjk-<l>j+i t k > (D 

Then an NxN square array of phase points has 2N(N- 1) 
connecting phase differences. 

For discrete subapertures, the errors on the phase 
estimates due to sensor noise or photon noise is well 
described as uncorrelated between subapertures. Thus 
the noise on the phase difference measurement S) k will 
be denoted n) k and will be assumed to have statistics 12 

^o^S^e,^ . (3) 

The actual measurements emerging from the wave- 
front sensor will be denoted by S l Jk and are the phase 
differences plus the noise error: 



(4) 



III. RECONSTRUCTOR MODEL 



The reconstructor is assumed to take the phase-dif- 
ference estimates S l jk and combine them linearly to es- 
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timate the wave-front phases <p Jk . Thus the most gen- 
eral form for the reconstructor algorithm is 



(5) 



where the * differentiates this b from a later symbol. 
The error in this estimator is 

^LXn + E b% l m „n l mn . 

\ tmn / Imn 

(6) 

We want to consider only estimators which give a 
perfect reconstruction when there is no noise (instan- 
taneously unbiased estimators). Thus one condition on 
the bfJn is 



(7) 



When there is measurement error, then the b% mn 
should be chosen to minimize the error effects in the 
# Jk . Here the choice is to minimize the mean- square 
error in 4> Jtt9 i. e. , to minimize 

(E«i)-n fl E (*&„) 2 , (8) 
where ( > indicates the statistical expectation value. 

IV. MINIMIZATION 

The general task is now to minimize where 



< 2 =(E.4,)=«oZ(&Z n ) 2 

\Jk / Imn 



Jk 



subject to the constraint that 



E &J*mn Smn 
Imn 



(9) 



(10) 



for any set of <f> Jk . 

For large arrays we assume that edge effects are 
small and that b*^ can be written as a function of j - m 
and & - n (translation invariance of &**£,„). Thus the 
solution to be derived applies to large arrays and will 
be perhaps suboptimal for finite arrays. It has also 
been demonstrated that the solution to be derived in the 
large array limit applies to a 2x2 array (small array 
limit). Thus both limits are correct and the deviation 
for intermediate arrays is likely small if nonzero. 
However, the question of optimality for a finite array 
remains open. 

Under this large array assumption, Eq„ (10) for the 
unbiased estimator becomes 



#i* =: E 6 /-m,*-n S mn > 



Imn 



(U) 



which is a convolution. 



Now we go to Fourier transform space since there 
convolution becomes simple multiplication. A A will de- 
note a Fourier transform: 



<?ifc = E&jfeSj* 9 



(12) 



where the different Fourier functions are defined by 



lm 
lm 

* Jk~ L*t ^ lm e 



(13) 
(14) 



(15) 



The minimization condition in Eq. (9) becomes (in 
Fourier space) minimizing in the following quantity: 



jkn 



(16) 



Now S) k is the difference between two <p m „ values de- 
fined by 

S/* = *rt-*i.w (18) 
In terms of Fourier components this is equivalent to 

§f k =(l-J>U"»)4> Jk . (20) 

Substituting these expressions for §) k in Eq. (12) 
gives a simple form of the unbiased condition in Fourier 
space : 

♦»-[S},(i-«* li/ *)+»J»(i-«* tt/ *)]i J , • (2D 

Since this equation is to be satisfied for any cp jk dis- 
tribution, the condition on b) k is 

Sj»(l-** ,i/ ^ + S} Jk (l-«* 4 * /J, )-l . (22) 

The minimization condition in Eq. (16) does not cou- 
ple different Fourier components, so it implies for b] k 
that the following expression is minimized for each (jk) 
value subject to the constraint of Eq. (22): 



6?» = (S}»)M$J») 1 . 



(23) 



The problem is now a simple application of variation- 
al calculus with the result that the optimum choice for 
b'n is 



(1 



~Zt{J/N 



kjk - 4 _ e ^U/ N _ e -2rij/ tf _ ^rik/ N _ e -2rik/ N > 



(24) 



(25) 



V. PHYSICAL INTERPRETATION 

This result for the optimum b l jh has a very simple in- 
terpretation when transformed back from Fourier space. 
3h Fourier space the estimated value of <j> Jk is found by 
substituting Eqs. (24) and (25) into Eq. (12) to get 

* >* " 4 Z e * 1 J/ " - e Zt ii/N - e 2 * ik/ N - e 2 * ihl N ' 1 ' 

Now multiply this equation by the denominator on the 
right side and realize that multiplying by an exponential 
in Fourier space simply translates the function in non- 
Fourier space. The result is 

$ jk [4 - e * u/N - €* ul *- e 2 "*'*- e* ikf *] 

= S) h (l - e* w N ) + 3}»<1 ~ e-* ik/ *) , (27) 
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N 


ib jkmn 


3 


0.683 


5 


0.726 


7 


0.754 


11 


0.807 


13 


0.825 


20 


0.870 



FIG. 1. Noise coefficient for 
various square arrays of size 
NxN. 



or in non-Fourier space this becomes 

Solving for <f> ife gives 

+ *(^-^-W + 55»-5f.«). (29) 
Now S) ft is the measured phase difference connecting 
<Pjk <l>j+i 9 ki a™* ^ft connects <p jk and <£ i)fe+1 [see Eqs. 
(l)-(2)]. Thus Eq. (29) has the interpretation that (p Jk 
is the average of its four nearest neighbors plus the 
average of the four connecting phase differences (with 
appropriate signs in the phase differences). 

VI. RECURSIVE IMPLEMENTATION 

The simple interpretation of Eq. (29) allows for a 
recursive implementation of Eq„ (29) where the new 
value <pj k of (p Jk is defined 



(30) 



The immediate questions are uniqueness and conver- 
gence of this recursive algorithm, and they can be 
demonstrated using standard linear algorithm tech- 
niques . 

VII. NOISE PROPAGATION 

Now that the convergence, existence, and uniqueness 
of a solution have been proven for the recursive algo- 
rithm, we return to the noise analysis. The mean- 
square error € 2 on the reconstruction phases was [from 
Eq, (9)] 



€ 2 



(31) 



jklmn 



where n Q was the mean-square error on a single phase 
difference measurement. The b] kmn have been calcu- 
lated using the recursive algorithm and the Yfi)^ is 
tabulated for various square arrays in Fig. 1. 

Note that for the arrays considered, the error prop- 
agation is less than unity, i. e. , less mean- square error 
appears on the reconstructed phases than goes in on the 
phase-difference measurements. A good fit to the data 
on the noise coefficient is 



^2 b)l mn ** 0.561 + 0. 103 IniST, 

jklmn 



(32) 



In addition, a simultation was constructed where ran- 
dom phase differences were generated and then used 
with the least-squares algorithm. The agreement on 
error propagation was within 5% of the analytical results 
in the two cases tested of N= 13 and N=2Q a 

VIII. NOISE CORRELATIONS 

Another aspect of noise error is its correlation 
across the aperture. This is significant for calculating 
optical transfer functions and optimal control algorithms 
for CI systems,, Here we calculate that correlation 
function 0 



The error on the phase estimate is 

Imn Imn 

which gives a correlation function of 



(33) 

(34) 
(35) 



Note that for r = s = 0 this becomes the noise error 
variance calculated in Sec« in. 

The b) k have been calculated for a 14x14 array of 
phase points using the recursive algorithm and the 
noise correlation function C rs plotted. It is very nearly 
rotationally symmetric so it is plotted in Fig. 2 as a 
function of radial spacing (r 2 +s 2 ) 1/2 . The correlation 
is seen to be quite long range, extending over roughly 
one -half the aperture. The correlation function is ap- 
proximately expressed 

C rs =C QQ e'* irZ + s2 > U2/1 * . (36) 

IX. HARDWARE IMPLEMENTATION 

The least-squares recursive algorithm has been im- 
plemented in digital and analog forms for compensated 
imaging systems. The algorithm requires a resistor 
net and N 2, op amps for analogue implementation and 
N z 16 -bit adders for digital operation. Note that divid- 
ing by 4 is a bit shift not a full division,, The errors 



Q, 




where N z = number of phase points in the square array 
(see Fig. 1). 



-0.2 



POINT SEPARATION V ? + s 2 
FIG. 2. Noise correlation function on the reconstructed 
phases. 
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induced by component tolerance and bit limitations have 
been seen to be negligible, and convergences times of 
around 10 us are reasonable for special-purpose, low- 
power hardware. 

X. CONCLUSION 

The optimum linear algorithm has been derived which 
converts noisy phase differences to absolute phase es- 
timates for a wave-front sensor* The algorithm has 
been given a simple, recursive form, and convergence 
and uniqueness have been shown. The noise propagation 
for this algorithm depends weakly on the size of the 
phase array (logarithmically) and is less than unity for 
the cases studied. This means that there is less rms 
error on the reconstructed phases than on the phase- 
difference measurements from the sensor. This error 
propagation is in marked contrast to a simple summing 
of phase differences which can produce many times 
more error. 

Hardware implementation, either in analog or digital 
form, is easily achieved and reconstruction times of 10 
lis are reasonable for low-power components. 
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Optimal wave-front estimation 

Richard H. Hudgin 
Itek Corporation, Lexington, Massachusetts 02173 

In a compensated imaging (CI) system the wave-front distortion is changing randomly in time and the data 
processing must follow the changes and even, to some extent, predict them. Furthermore, the system is not 
infinitely fast; it has time delays and bandwidth limitations which can be quite significant. Couple to this the 
substantial photon noise in the wave-front measurements, and it becomes desirable to find the optimal data 
processing to use for system control. This paper derives the optimal linear control algorithm assuming the 
aberration and noise statistics are known for a compensated imaging system, and discusses the properties of 
this algorithm. 



I. INTRODUCTION 

In a compensated imaging (CI) system the wave-front 
distortion is changing randomly in time and the data 
processing must follow the changes and even, to some 
extent, predict them. Furthermore, the system is not 
infinitely fast; it has time delays and bandwidth limita- 
tions which can be quite significant. Couple to this the 
substantial photon noise in the wave-front measure- 
ments, and it becomes desirable to find the optimal 
data processing to use for system control. This paper 
derives the optimal linear control algorithm, assuming 
the aberration and noise statistics are known for a com- 
pensated imagining system, and discusses the prop- 
erties of this algorithm . 



II. SYSTEM MODEL 

A d system consists of three main components: a 
wave-front sensor, a data processor, and an active 
mirror. The wave-front sensor produces photon noisy 
wave-front estimates at every time interval 6 and in- 
puts them into the data processor . The processor, in 
turn, combines the measurements linearly and gene- 
rates a predicted wave front which is sent as com- 
mands to the active mirror. 

Thus we will assume a set of N phase estimation 
points at positions Xj:j=l~N 9 and the wave-front sen- 
sor will generate noisy phase measurements 4> q ^l J9 
t-Kd) for K=Kq, -oo. 
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The value of K Q represents the time delay of the sys- 
tem, measured in units of the sensor interval 6. An 
ultimate system would have K 0 = 1, (i.e., the data is 
used right after it is measured with a mean delay of 
one cycle), but the data processing requirements may 
be too heavy for that limit to apply. Thus iT 0 = 2 or 3 
may be more meaningful for an operational system, and 
the performance of the estimator as a function of the 
delay K Q 8 is studied below. For a postprocessing esti- 
mate of the wave front, all data can be used, so K Q 
= - oo in that case. Actually, there will be a finite cut- 
off in the number of wave fronts processed due to CPU 
limitations even for postprocessing estimates of wave- 
front error. The effects of restricting the number of 
terms are also studied. 

The dataprocessor takes these noisy wave-front 
estimates 0 O (X/, t-K8) and combines them linearly 
with weights A {jk to get an estimate $ (x f , t) for the 
wave-front at time t: 



$(x i9 t)^J2A is J 0 (x j9 t-kd) . 



(1) 



Note that the A iJk will usually be different for each x. . 

The error of this estimate is the difference between 
the estimate $(x t , *)and the actual distortion <f> a (x it t). 
The mean-square error E will be calculated below and 
depends on the weights A lJk , the wave-front statistics, 
the cycle time 5, the time delay K 0 5, and the sensor 
noise w(x / , t- kd) in the measurements $ Q (x s , t- kd). 

Perhaps we might mention here that several sources 
of error are not considered in this analysis. The error 
due to discretely sampling the wave front is not cal- 
calculated, and is treated elsewhere. 12 Also, errors due 
to atmospheric dispersion, 13 mis ceUaneous system opti- 
cal distortions, and isoplanatism are not considered. 
The final system error will be a sum of these different 
error terms, so the present paper is by no means a 
full system analysis . 

III. MEAN-SQUARE ERROR 

The error variance E of the estimate $ can be cal- 
culated from the statistics of the wave-front error and 
the sensor noise. We choose to characterize the wave- 
front statistics by a structure function D defined by 

*) = <[<Uy+x, f + *)-0 a (y, t')] z ), (2) 

where ( } denotes statistical expectation value. 

The sensor noise is assumed to have a known cor- 
relation function C jklm defined by 

Cwm=(n(Xj, t-k6)n(Xi, f-m6)> . (3) 

In the case of a wave-front sensor which is a shearing 
interferometer or a Hartmann sensor, the function 
c mm has been analyzed 4 and is approximately expressed 
as 

where D 0 is the diameter of the telescope aperture and 
C 0 is the mean-square noise. Note that measurements 
at different times (i.e., K±m) have independent pho- 



tons and thus independent noise errors. 

Using these statistics in Eqs. (2) and (3), one can 
calculate the error variance E as follows: 

*i»<[*.fc<, *)]"> (5) 

= j£ AiMxj,t-k6)J^ , (6) 

where we substituted Eq. (1) for <£. 

Now express <£ 0 as <p a plus the noise error n to get 



f n(x 



(7) 

It is actually not the absolute phase that is mea- 
surable in the aperture, but relative phase. Adding a 
uniform phase to the aperture is not physically mean- 
ingful- thus the mean-square error E { should involve 
wave fronts with zero mean, and to get the physically 
meaningful error, we interpret <f> a as a wave front of 
zero mean. Subtracting this mean allows the use of the 
wave-front structure function in calculating £ r and 
gives the following result; 

- Z f-*0) + n(x Jf f-JW)] 

■ j7 t Yj^^uI^Mx, t-kd) + n(x t , f-fe6)]J^ 



fc=k n ~<x> 



-- E i*+ Z >A iJk B Jk+ £ A {Jk A itm C mm , 
/si-* /, 1*1- 

*=*0-« ft, m=fc 0 -» 



(8) 
(9) 



where 



W 2 



Z {D(Xj-*-i, 0)+C mo }, 



u ux-n 



(10) 



£ M = D(x,-x y , kd)+C {kJQ 

-jjf H {D(X|-x 7 , k5)+D(x J -x l , -fe6) 

1 V 

+ c mo+C iOU }+— l 22 {D(x l -x m ,k5)+C t > m0 }, 
Cmm =-iD[x,-x l , (m- k)6] - I C mn 

+ 2F E {Dhi-X,, (m-A)5] + i?[ Xi - Xn , 

(m-k)5]+C„ klm+ C Jlmm } 
-2&* 5 {D[x n -x„ (w-^Sj + C^J. (12) 



(ID 



Now Eq. (9) is aquadratic equation involving the A iJk . 
The values ofA^ which minimize^ can be found by setting 
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the variational derivative of E { with respect to A ijh equal 
to zero. This produces the following set of linear, 
simultaneous equations which can be solved for the 



lm 

IV. TEMPORAL AVERAGING 



(13) 



The first attempt to use these equations was to take 
the past wave fronts at a single point and to estimate the 
correction from them. Thus only temporal weighting, 
not spatial weighting is involved, which implies 



• 0 for j±i 



(14) 



Here we set K 0 - 2 to represent a one cycle process- 
ing delay in the system and found the error as a func- 
tion of sensor noise and 5 for one, two, and three 
terms in the past (with optimum weights). Including 
more than three terms caused little improvement in 
the range studied. The results are presented in nor- 
malized form in Fig. 1. 

We have assumed a Kolmogorov atmospheric structure 
function where 



D(0 9 t) = d Q t 5f3 (waves 2 ), 
and have used the noise correlation function 
C /Wnt =C 0 5 feM (waves 2 ) , 



(15) 



(16) 



where each wave-front measurement has its own inde- 
pendent noise error. 

If there is no photon noise, then the error using a 
single term which is two cycles old is d 0 (25) 5/3 . This 
is the y intercept of the "one wave-front" curve. Using 
more terms is always better, but only slightly until the 
photon noise gets very large compared to the time de- 
lay error for one cycle, i.e., 

C 0 »rf 0 6 5 ' 3 . (17) 

When the photon noise is that large, then a larger 6 
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FIG, 1. Dependence of the estimation error on sensor noise 
and number of terms for temporal averaging. 
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FIG. 2. Atmospheric turbulence distribution used to evaluate 
the estimator gives r Q = !0 c (equal turbulence in tropopause 
peak and the exponential) . 



can be used without significant degradation. Thus a two- 
term predictor is nearly optimal when spatial correla- 
tions are ignored. 

V. SPACE-TIME WEIGHTING 

The residual error of the optimum space-time esti- 
mator was studied next. We chose an array of phase 
points that was 21X21 and used 40 wave fronts in the 
past. The 200 phase points with the largest correla- 
tion to the point to be estimated are chosen out of the 
2ix 21 X40= 18 081 candidates. Then the 200 optimal 
weights are calculated. The resulting mean -square 
error is felt to be a close approximation to the mini- 
mum possible error. 

When the effects of fewer terms were explored, the 
N terms with the largest weights among the 200 are 
used. The optimum weights are recalculated for these 
N terms, and then the mean- square error is derived. 

The dependence of the error on position in the 21 
x21 aperture was also studied as well as the selection 
pattern of the terms. 

The wave-front statistics were derived from a 20- 
layer atmosphere with an exponentially decaying low- 
level turbulence and a strong tropopause layer of tur- 
bulence. The turbulence strength C\ is plotted as a 
function of height in Fig. 2, and corresponds to a co- 
herence length of r 0 = 10 cm (i. e. , 1 arc sec seeing). 
The structure function used was 

*> = "\r^ j dh(^ n (h)\x~hQ}t\^ 3 (waves 2 ), (18) 

where cj is the angular slew velocity and X is the mean 
wavelength of light. This corresponds to a Kolmogorov 
type of turbulence which moves due to the slewing of 
the line of sight through the atmosphere. 

The noise correlation function used was the one in 
Eq. (4). 

VI. OPTIMUM SPACE-TIME ESTIMATION 

The results of the 200-term estimator are shown in 
Fig. 3 for various levels of sensor noise and system 
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WAVE-FRONT 
NUMBER 



2 4 6 8 10 12 

RELATIVE SENSOR NOISE {C 0 /d 0 5 

FIG. 3. Dependence of estimation error (100 terms) on system 
speed and sensor noise. 

delay. The slew rate corresponds to the rotational 
velocity of the earth, and the sensor interval is 5 = 35 
ms. The error increases as sensor noise increases 
and for a real-time system does not go to zero (even 
without sensor noise) due to a time delay of at least 
one- cycle 5 between when the data is measured and 
when it is used. This delay error is gone as soon as 
the postprocessing region is entered {K Q < 0). There 
the wave front to be estimated is observed, and thus 
the estimation error vanishes when sensor noise 
vanishes. 

If there were no sensor noise, the error for a practi- 
cal system with i£" 0 =2 (i.e., one cycle interval 6 for 
data processing) would be d 0 (26) s/z if only the last 
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FIG. 4. Predictor error variance vs relative sensor noise and 
number of terms QT 0 = 2). 
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FIG. 5. Term ordering by weight for the optimal estimator. 

measurement were used in the estimator. The reduc- 
tion in error by a factor of 10 from the one -term esti- 
mator to the optimal estimator indicates the importance 
of using more than one term . 

A posteriori estimates (here approximated by K Q 
= - 19 which used 19 past and 20 future wave fronts to 
estimate a current wave front) are substantially better 
than the real-time estimates. They can be used for 
possible postprocessing applications 

VII. DEPENDENCE ON THE NUMBER OF TERMS 

In Fig. 3 the contrast for jRC 0 = 2 between one term and 
200 terms was considerable. In Fig. 4, the depen- 
dence of the estimator error on the number of terms 
for a practical system {K Q = 2) is shown. For sensor 
noise smaller than or equal to the error due to the 
system time delays, roughly 20 terms will give nearly 
optimum performance. As the sensor noise increases, 
more averaging is required so that more terms are 
useful. If sensor noise increased beyond twice the 
system delay error, then using a longer sensor inte- 
gration time is probably preferable to using more 
terms. 

VIII. TERM SELECTION 

In Fig. 5 the order of the first 30 terms is tabulated 
by the magnitudes of their weights. The selection pat- 
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TURBULENCE FLOW 




and applied to various cases. In general, a substantial 
improvement results from using more than the last 
sensor measurement. Spatial correlations in the wave 
front are seen to be significant for Kolmogorov turbu- 
lence and the estimation is very poor v:hen they are 
ignored. 

For finite apertures, the accuracy of the estimator 
is nearly independent of position except at the leading 
edge of the aperture where the new turbulence is ap- 
pearing. It may be desirable to mask off the leading 
edge in a practical CI system. 




FIG. 6. Estimation error as a function of position in the aper- 
ture. 

tern follows the flow the turbulence back in time and 
spreads out only very slightly orthogonal to the turbu- 
lence flow. The first wave front available is two be- 
hind the one being estimated for a K Q = 2 system. 

IX. SPATIAL DEPENDENCE OF THE ESTIMATOR 

One question to consider in designing a CI system is 
how the error varies across the aperture. If the edges 
of the wave front are significantly worse than the in- 
terior, then masking off the aperture before imagining 
may improve the system performance. 

The numerical analysis has been performed for a 
variety of locations in an aperture with a 21 x 21 point 
grid with the results presented in Fig. 6. There is a 
slight increase in estimation error near the edges of 
the aperture, but it is only significant right at the lead- 
ing edge when the turbulence is just entering. Conse- 
quently, only the leading edge would be a candidate for 
masking off. The results are presented for two levels 
of sensor noise with the same general result. 

X. CONCLUSION 

The optimum linear estimator for a compensated 
imagining system with sensor noise has been developed 



The number of terms required to get nearly optimal 
performance in the estimator increases as the sensor 
noise increases, but for sensor noise less than or equal 
to the error caused by system time delay, 20 terms is 
nearly optimal. 

Postprocessing estimates of the wave front are sub- 
stantially more accurate than real-time estimates 
which may be useful in image enhancement. 
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Phase estimates from slope-type wave-front sensors 
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In a recent issue of this journal Fried and Hudgin 1 ' 2 show 
that the expected mean-square wave-front error due to mea- 
surement noise over a square array of phase-difference-type 
wave-front measurements depends on the total number of 
measurement points N according to A + B IniV. The methods 
employed involve a computer least-square solution of the 
equation net. In this Letter, it is shown that the ln/V depen- 
dence can be derived analytically, and suggested that the 
least-square solutions are similar to solutions of Poisson's 
equation with Neumann boundary conditions. 

Minimizing the mean-square wave-front error for a com- 
pensated imaging system with a phase difference wave -front 
sensor has been shown 3 to be consistent with a wave-front 
corrector that solves Poisson's equation. Since only phase 
slopes are measured, the boundary condition involves only the 
normal derivative (Neumann boundary conditions). Solving 
Poisson's equation is a classical problem which has a wide 
variety of solution techniques as discussed by Morse and 
Feshbach, 4 among them finite difference solution methods. 
With Neumann boundary conditions the phase can be ob- 
tained to within an arbitrary phase provided that the 
boundary slopes satisfy an integral requirement like 

/s(r) . n dl = 0, (1) 

where f denotes integration around the boundary, n is a 
normal to the boundary, and s is the slope at point r. A 
knowledge of the Green function for the measurement con- 
figuration is all that is required to find a solution. On a lattice 
denoted by (m,l) the phase $ can be obtained by 

Q,t 

+ Z(sB*n)g(qB,tB\™>D> (2) 

B 

where s# • n is the normal component of the slope measure- 
ment at the boundary and g is a Green function that has a zero 
normal derivative at the boundary. The quantity F is simply 



the result of taking the divergence of the slopes, 

V - s = -F. (3) 

The solution to the phase-difference problem therefore cen- 
ters solely on finding the appropriate Green function g. 

A simple analytic Green function exists for the infinite 
boundary problem. Let us consider this Green function and 
neglect the boundary term. For convenience let us use a 
continuous solution rather than the discrete form given by Eq. 
(2); that is, 

*(r) = - Jdr'[V-s(r')k(r|rO. (4) 

Taking the Fourier transform of Eq. (4) and noting that the 
infinite boundary Green function transform takes the form 

G(k|k0 = (l/47r 2 fe 2 )6(k-k0 (5) 

leads to the Fourier equivalent of Eq. (4), 

*(k) = -ik.S(k)/27rfc 2 . (6) 

This form can be used for calculating the mean-square 
wave-front error due to noise. Multiplying Eq. (6) by <£* (k ) 
and taking an ensemble average gives 

<4>*(k)*(k)) = {[k.S*(k)][k.S(k)])/(27r/2 2 ) 2 . (7) 

Now make the usual assumption that the slope noise due to 
the measurement sensor is white over the system pupil and 
that the x-y components are independent so that the spatial 
power spectral density N 0 of the wave-front slope measure- 
ments is defined by 

N 0 = %<S*(k)-S(k)>. (8) 

The spatial power spectral density of the wave front given by 
Eq. (7) is thus 
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<**(k)*(k)> =AT 0 /2ir 2 * 2 . (9) 

The singularity at the origin in Eq. (9) is a result of the fact 
that a solution can be obtained to within an arbitrary phase. 
Taking the constant phase over the pupil to be zero is equiv- 
alent to multiplying Eq. (9) by 1 - Af(fe), where 

A?(ft) = sinc 2 (7r/z x D) sinc 2 (7r/e y D) 

(square pupil of size D ) , (10) 

A?(A) = AJ\(2irkR)y{2irkR)^ 

(circular pupil of radius R). (11) 

Thus the mean-square wave-front error due to measurement 
slope noise E N can be written 

EN .!kj (12) 

2-lZ 1 jNyquist ft 2 

where the limit on the integration is the spatial Nyquist fre- 
quency K, which is related to the total number of sampling 
points, N, in the pupil of dimension D by K = VN/2D. 
Numerical integration of Eq. (12) yields 

Square aperture: E N = 2N 0 [0.0536 + 0.0795 IniV], (13) 

Circular aperture: En = 2N 0 [0.0068 + 0.0796 lnN]. (14) 

Equation (13) should be compared with Fried's result, which 



in the present notation is 

En = 2N 0 [0.3279 + 0.0815 lniV]. (15) 

The only significant discrepancy between Eqs. (13) and (15) 
is in the constant term which we take to be due to our ne- 
glecting the boundary term in Eq. (2). It should be noted that 
aberration errors like tilt and astigmatism contribute only 
through the boundary term. This is because V • s = 0 for tilt 
and astigmatism inside the pupil. Thus neglecting the 
boundary term results in omission of significant low-order 
aberration errors. The origin of the IniV dependence is a 
consequence of the Ilk 2 power density of the wave-front 
phase, i.e., of the two-dimensional integration of the wave- 
front slopes. Put another way, white wave-front slope noise 
produces Ilk 2 wave-fronts errors, which means larger errors 
associated with low spatial-frequency distortions. 

I believe that the relaxation method used by Fried for a 
square lattice to obtain the mean-square wave front is similar 
to a method suggested by Morse and Feshbach 4 for finding 
the Green function. Because Eq. (2) applies for arbitrarily 
shaped pupils, the method should be of general usefulness in 
obtaining the phase from slope measurements. Ail that is 
required is the determination of the Green function corre- 
sponding to the specific geometry of the problem. 

1 D. L. Fried, J. Opt. Soc. Am. 67, 370 (1977). 
2 R. H. Hudgin, J. Opt. Soc. Am. 67, 375(1977). 
3 R. J. Noll, Proc. SPIE 75, 39 (1976). 

4 P. M. Morse and H. Feshbach, Methods of Theoretical Physics, I 
(McGraw-Hill, New York, 1953), p. 698. 



Brightness exponent for the periphery in the Bloch region 

Naoyuki Osaka 

Department of Psychology, Otemon-Gakuin University, Ibaraki, Osaka 567, Japan 
(Received 11 August 1977) 

Using a method of direct magnitude estimation, perceived brightness was measured in the dark-adapt- 
ed eye with brief flashes of varying duration in the Bloch region (1-100 ms) and retinal lock (0°-40°) for 
the pho topic luminance levels covering the range between 140 and 0.14 cd/m 2 in steps of 1 log unit. Per- 
ceived brightness increased as a function of flash luminous energy (product of luminance and duration) 
up to critical duration of approximately 100 ms. The brightness power exponent for brief flash was 
found to be 0.48 ± 0,01 in the fovea, whereas about 0.44 ± 0.01 in the periphery. 



It has been established in threshold studies as Bloch's 
law that luminance power is integrated over time up to a 
critical flash duration (t c ) of approximately 100 ms (Bloch 
region), beyond which temporal summation ceases and the 
threshold is then defined solely in terms of flash luminance 
(L). Bloch's reciprocity law (Lt = const.) and a decrease in 
t c with increasing L, in the fovea, have also been established 
for suprathreshold moderate luminance range in the dark- 
adapted eye. Perceived brightness grows as luminous-energy 
(Lt ) raised to a power in the Bloch region. Brightness can be 
expressed as power law of the form \p = k 2J0 ±S - E -, where \p 9 E, 
j8 ± S.E., and k indicate perceived brightness, luminous en- 



ergy, exponent ± standard error, and a scale factor, respec- 
tively. 

In the foueal viewing, the brightness exponent has been 
found to be approximately 0.50 for brief flashes below the t c , 
since the exponent has been estimated approximately 1% 
times larger for transient-state flashes than for steady-state 
flashes. 1 " 3 If the exponent for longer flashes (out of the Bloch 
region) is assumed to be 0.33, 1 * 3 the exponent for brief flashes 
(within the Bloch region) may be estimated to be close to 
0.50. 2,3 The study of Raab 4 gave a flash exponent of 0.5, that 
of Mansfield, 2 0.50 ± 0.03, that of Stevens and Hall,* 5 0.45, and 
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Comparison of wavefront sensor configurations 
using optimal reconstruction and correction 

Edward P. Wallner 

Itek Optical Systems, a Division of Itek Corporation 
10 Maguire Road, Lexington, Massachusetts 02173 

Abstract 

Recent work on optimal reconstruction of wavefronts from slope measurements has developed 
methods based on more realistic models of sensor devices. These methods also resolve cer- 
tain ambiguities present in earlier reconstruction methods. In the present paper , these 
methods are applied to sensor configurations which have been proposed in the literature and 
the errors in estimating and correcting wavefronts are compared. 

Current wavefront reconstruction methods 

Most means of defining the shape of a wavefront are based on measurements of local wave- 
front slope. For instance, in a shearing interferometer, two displaced copies of the wave- 
front are interfered and the phase of the interference pattern is a measure of wavefront 
slope in the direction of the displacement or shear. In a Hartmann test, the displacement 
of a spot from its nominal position is also a measure of local slope. 

The wavefront itself must be reconstructed from the slope data to evaluate or correct its 
errors. One scheme for this was described by Rimmer. 1 Here the measurements are made by 
arrays of x and y slope sensors which are displaced as shown in Figure 1. The output of 
each sensor then approximates the phase difference between the grid points on the edges of 
the sensor. This ties the phases at the grid points together in a network such as that in 
Figure 2, which corresponds to the configuration of Figure 1. There are many redundant 
paths in this network allowing the measurement noise to i>e reduced by averaging. 
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Figure 1. Displaced subaperture sensor 
configuration. 



Figure 2. Network for phase reconstruction 
with displaced subapertures. 



396 



If only measurement noise is considered and it is equal and independent in each sensor, 
the phase estimate giving the least mean square error is the average of the estimates 
obtained from the phase at adjacent points plus the measured phase differences. 1 ' 2 / 3 



(1) 



where 



4>j = phase at 



j-th grid point 



^ij ~ measured phase difference between grid points i and j 



The solution to this system of equations is generally found by iteration, and it may be 
expressed as a matrix that multiplies the vector of phase differences to produce the 
vector of phases 



>j = M jk 



(2) 



where 



M jk = solution matrix 



A<|> k = k-th phase difference measurement 



The average phase over the aperture is set to zero, since a constant phase will neither 
affect optical performance nor be detected by the slope sensor. 

A different wavefront sensor configuration, representing a Hartmann type of sensor, was 
analyzed by Fried. 4 This sensor, as shown in Figure 3, measures the x and y slopes in the 
same subaperture. If the measurements are resolved onto axes rotated by 45° (or if the 
measurements are made in those coordinates) , the new measurements approximate the phase 
differences between the corners of the subapertures. Figure 4 shows the resulting network 
for the configuration in Figure 3. There are now two interlaced grids of points that 
share no common grid point. The phase can be reconstructed for each net as before, but 
each has an arbitrary constant phase term. The overall average phase can again be set to 
zero but the difference between the two nets remains as a degree of freedom. 5 A difference 
in the average phases will lead to a "checker board" phase pattern to which the sensor is 
insensitive . 
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Figure 3. Common subaperture sensor 
configuration . 



Figure 4. Network for phase reconstruction 
with common subapertures. 
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In addition to the arbitrary degree of freedom in the latter case, there are several 
other limitations to this approach to wavefront reconstruction: 

1. The slope measurement is treated as a point-to-point phase difference rather than 
an average over a finite area of the aperture; . 

2. Estimates are obtained for the grid points only with no guidance on how to interpo 
late between them; and 

3. No use is made of the statistics of the wavefront being estimated. 

In the next section, a reconstruction method that overcomes these limitations is 
developed. 

Optimal wavefront reconstruction from slope measurements 

The system considered here consists of an aperture and a number of wavefront sensors, 
each of which measures the weighted average of wavefront slope over some region within the 
aperture. The aim is to use the measurements produced by these sensors to construct an 
estimate of wavefront phase throughout the aperture with minimum mean squared error. 

The aperture of the optical system can be described by a weighting function proportional 
to the intensity of light at the wavefront. It is convenient to normalize this function 
so that 



/: 



dx W.(x) = 1 



(3) 



where x = two-dimensional vector position in the aperture plane 

W A (x) = aperture weighting function 

f dx = indicates integration over the entire aperture plane 

J — CD 

The wavefront phase or optical path difference is described by the function ifj (x) . Because 
any uniform phase over the aperture does not affect optical performance, we can equiva- 
lently deal with the phase with aperture average removed, <j> (x) . For this function 



j dx W A (x) <Mx) 



(4) 



where * (x) -♦(x) - j dx'W A (x 1 ) Y(x') 



i|Mx) = input wavefront phase 



Each wavefront sensor forms a weighted average of the wavefront slope over some region of 
the aperture. This average is corrupted by additive noise. The signal out of the n-th 
sensor is then 



dx W sn (x) [ $ n (x) + v(x)l (5) 



where s R = signal from n-th sensor 

W (x) ■ weighting function for n-th sensor 
sn 

<f> n (x) ■ slope of wavefront in the direction of sensitivity of n-th sensor 

v(x) = noise distribution over aperture 

The signal can be expressed in terms of the wavefront phase by integrating Eq. (5) by 
parts : 



s n = f° dx r-Wg n (x) 4>(x) + W sn (x) v(x)I 



(6) 
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In Eq. (6) use has been made of the fact that W sn (x) goes to zero at infinity. The noise is 
assumed to be unbiased, and uncorrelated with wavefront phase 

<v(x)> = <v(x) <Mx')> = 0 (7) 

The estimate of phase at some point x in the aperture is taken as a linear combination of 
the slope measurements 

<f>(x) = £ s n r n (x) (8) 
n 

where (j> (x) = estimated phase at x 

^ n (x) = n-th weighting function in estimate 

We wish to determine a set of weighting functions r (x) that minimize the error in the 
estimate. The error is 

e(x) = (j) (x) - 4>(x) 

= E s n r n (x) - <Mx) (9) 
n 

The mean square error in the estimate is then 

<e 2 (i)> = Z l r (x) r . (x ) <ss 

~ « - n n n n 

n n 

- E r n (x) <s n (j)(5)> 
n 

- r n -*(x) <s -*<^) > + <4> 2 (x)> (10) 
n " 



where <f> = indicates the expected value of f 

Eq. (10) can be simplified by defining 

S I = <s s ^> 
nn n n 



f ax- f ax"[w; (x-j w; n .(x") «Mx*) *(x-)> 

•'-co — CO [__ 

+ W sn (x') W sn . (x") <v(?) v(x")>] (11) 



f n (x) = <s n 4)(x)> 



- J dx' W^ n (x') «f>(x) *(x')> (12) 

Eq. (7) has been used in reducing Eqs. (11) and (12) . 

Substituting into Eq. (10) gives the simplified expression for the mean square error at x 

<e 2 (x)> « I £ r (x) r . (x) S,. 

_ _ ■* n n nn 

n n 

- I r B (5) f n (x) 

- r n .(5) f n -(x) + «j> 2 (x)> (13) 
n ** 
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The optimum value of r n (x) is found by differentiating Eq. (13) with respect to r^x) and 
equating to zero: 



d <e^(x)>/dr .(x) - I r n (x) S , - V (x) 

n 

= 0 



(14) 



The solution for the optimum estimation function is found by inverting Eq. (14). 
Indicating the optimum by a superscript * yields 

r*(x) = I S~K f n ^(x) (15) 

n nn n 
n 

where - optimum estimation function for n-th sensor 

S -1 ^ - inverse of , matrix 
nn nn 

The matrix S nn - will not be singular unless there are redundant slope measurements with the 
same noise , which would not occur in any real system. However, as the number of sensors is 
increased and the noise is reduced, the matrix becomes less well conditioned and special 
inversion techniques may be required to maintain accuracy in computation. 

The resulting minimum mean square error at x is found by substituting in Eq . (13) : 

<e 2 (x)>* = <<|> 2 (x)> - Z r*(x) f(x) (16) 

n 

Finally, the minimum mean square error averaged over the aperture is 

dx W„ (x) <e 2 (x) >* (17) 



•i-JL 



The accuracy of the estimate will depend on the signal-to-noise ratio of the measurements. 
With very noisy data, the matrix S nn - becomes large and the estimated wavefront becomes 
small as the a priori estimate of zero phase becomes more reliable than the measurements. 
In this limit, the error approaches the original wavefront and is independent of the density 
of sensors. 

For small measurement errors the matrix S nn ^ approaches a fixed value and the estimation 
error approaches a limit set by the inability of the finite array of sensors to measure or 
infer the wavefront perfectly. 

Representation by structure function 

In the case of atmospheric turbulence, the statistics of the input wavefront are defined 
in terms of a structure function 

D(x\ x') = <t*(x) - iMx")] 2 > (18) 

2 — 

The correlation of phases with averages removed required in the computation of <<f> (x)>, 
S and f (x) can be derived directly using Eq. (18). The average removed phase is 

nn n 

<Mx) = j dx" W A (x") [<M*) - i|>(x")] (19) 
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The correlation required is 
<<f>(x) <Mx")> 



= J a d *" J W A (x") w A (x-') <[<p(x) - t(x")] [t(x') - *(x"")]: 



ax'" w A (x") w A (x'") 



-00 ^_oc 

[D(x, x") - D(x, x"") - D(?', x-) + D(x", x'")] 
= D(x, x") + g(x) + g(x') - a, {20 ) 



where 



g(x) = h j dx" W A < X> ') D(x, x") 

dx"" J ax'" W A (x — ) W A (x"') D(x", -) 



(21) 



a = 



(x) g(x) (22) 



The mean square wavefront phase at x is 

<* 2 (x)> = 2g(x) - a 
The aperture average of the mean square phase is 



(23) 



/ ax w a 



(x) <<j> 2 (x) > = a 



(24) 



In the evaluation of S by substituting in Eq. (11) all terms on the right in Eq. (20) 
except the first can be seen to integrate to zero, leaving 

= f f dx- [-h W sn (5l W^U") D(?, x-) 



S nn- 



+ W sn (x') W sn .(x-) <v(x') v(x")>] ( 25) 
The function f n (x) is evaluated by substituting Eq . (20) in Eq . (12). 

f n (x) = j dSr W^ n (x-) [^D(x, x') - g(?)] 

= H £„ "m 5 ' 1 D( *' *' } " C n (26) 

The constant c n can be seen to be the value of the first term on the right averaged over 
the aperture, so that the function f n (x) has no net average phase. The same will also be 
true of the estimation functions r n (x) . 

Photon noise 

In many systems of interest, the limit on performance is set by the noise due to photon 
statistics in the wavefront measurements. The photon noise for a single wavefront can be 
represented by making v(x) a spatially white noise function. The correlation of the noise 
on two wavefront measurements depends in detail on the configuration of the wavefront 
sensor used. If the two measurements are made on separate wavefronts, either in sequence 
or separated by a beam splitter, they will be uncorrelated even though they refer to the 
same area in the aperture. 
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Measurements of orthogonal slopes made simultaneously on the same wavefront will also 
be uncorrelated. To reflect these conditions, the correlation of the noise functions may 
be written 

<v(x) v(?)> = k nn . 6(x - x") (27) 
2 

where a R = photon noise density 

k „ = T n • T n , for n and n' measured simultaneously on the same area of the same 
nn wavefront 

= 0 otherwise 

l n = unit vector in direction of slope sensitivity of n-th measurement 
6(x) = Dirac delta function 
The noise term in Eqs. (11) or (25) may then be written 



dx'" W (x-) W „(x") <v(x') v(x")> 
sn sn 



= k„ . al f dx' W (x') W en .(x') (28) 
nn n J_ m sn sn 

In most cases the basic measurements will be disjoint and 3c nn > would be the Kronecker delta 

if each measurement is treated separately. It may be convenient to combine a number of the 

basic measurements before processing in order to reduce the order of the matrices involved, 

and in this case the more general form of k „ would be needed. 

nn 

Optimal wavefront correction 

The analysis above applies directly to the problem of how best to correct a wavefront 

given the outputs of a set of wavefront sensors. A If each sensor drives an actuator that 

has a response of r* (x) , the total correction of i = I s r* (x) will give the best possible 
n n n 

correction for a linear system. No additional degrees of freedom beyond the number of 

sensors can improve the correction. 

The commands to the correction functions developed here are not in general statistically 
independent so that if any corrections are omitted, those remaining would not give an 
optimal correction with the reduced number of degrees of freedom. By properly transforming 
the commands, however they can be made independent, allowing the best correction for any 
number of degrees of freedom to be found. 

If G „ is an orthogonal matrix consisting of unit eigenvectors of the matrix, the 

transformed commands will be 

s n = J, G nn- s n' < 29 > 

and then optimum transformed correction functions will be 

r*'(x) « Z G r* (x) (30) 
n n ^ nn n 

The reduction in the mean square residual error averaged over the aperture attributable to 
the n-th transformed correction is 



A n a£ = J dx W A (x) f^(x) r*'(x) (31) 



where f*(x) = E G nn . f R (x) (32) 

n * 

3y adding correction functions in the order of magnitude of their overall error reduction, 
the minimum expected residual error is obtained for any number of degrees of freedom. 
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With low noise and well distributed sensors, the lower order correction functions will* 
be similar to those of the Karhunen-Loeve expansion that corresponds to the case of con- 
tinuous, noise-free measurements of the wavefront. This expansion has been developed for 
a circular aperture and a Kolmogorov structure function by Wang and Markey.6 

The case of the best possible correction given a specified set of non-optimal correctors 
has been treated elsewhere." 7 

Computational results 

The theory derived here has been applied to a unit square aperture with uniform illumi- 
nation with wavefront sensors that measure the average slope with uniform weight over square 
subapertures filling the aperture. The measurement noise was taken to be photon noise. 

Three configurations of. sensors were considered corresponding to the three cases treated 
by Southwell. 8 These are shown in Figure 5 for the case in which the subaperture dimension 
L is half the aperture dimension A. In the first configuration, as shown in Figures 1 and 
5a, x and y slopes were measured in displaced apertures. The sensors at an edge parallel 
to the direction of the slope measurement have half of their area outside the aperture, and 
the signal and photon noise power are accordingly halved for them. 



X SLOPE SENSORS 




(A) CONFIGURATION 1 




(B) CONFIGURATION 2 (C) CONFIGURATION 3 



Figure 5. Sensor configurations for L = A/2. 
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In the second configuration, as shown in Figures 3 and 5b , x and y slopes are measured 
in the same subaperture. Here all sensors have the same signal strength and noise power. 

The third configuration as shown in Figure 5c also measures x and y slopes in the same 
subapertures, but the array of sensors is shifted by half a subaperture in the x and y 
directions. Sensors at the edges of the aperture have their signal and noise powers 
halved, and for those in the corners of the aperture the signal and noise powers are 
quartered. 

The wavefront being corrected is described by a Kolmogorov structure function 
D(x) = 6.884 |x/r 0 | 5/3 (33) 
where D (x) = structure function in radian squared of phase 
r^ = Fried 1 s coherence length (m) . 

The mean square uncorrected phase error over the full aperture for this structure function 
is 



2 = 1.3103 (A/r n ) 5/3 (34) 
w u 

2 2 
where = phase error (rad ) 

A - dimension of side of aperture (m) 
The mean square slope measured over a sensor subaperture is 

a 2 = 6. 4051 L~ 1/3 r ~ 5/3 (rad 2 irf 2 ) (35) 
s u 

where L = side of sensor subaperture (m) 

The photon noise on the measurement expressed in radians squared of phase difference 
is inversely proportional to the number of photons used in the measurement 

a 2 = K/BL 2 (rad 2 ) (36) 

2 2 
where a = sensor noise (rad ) 
n 

2 -2 -1 

K = constant (rad photons m s ) 

-2 -1 

B = wavefront irradiance (photons m s ) 

The constant K will depend on the distribution of radiance over the source of the wavefront, 
the effective exposure time used in the measurement, the optical efficiency of the system, 
and the details of operation of the wavefront sensor. 

Figure 6 shows the shape of the correction functions and the mean square residual error 
for an array of sensors in the second configuration with L = A/2. In this case, the four 
x and the four y correction functions have the same shape except for reflections and rota- 
tions. The mean square residual error represents the result of applying the corrections 
to a series of random wavefronts, squaring the residual phases after correction, and 
averaging over the series. 

When the noise is very high, the correction function approaches a pure tilt, but the 
shape is unimportant because the magnitude of the correction becomes small. The residual 
error then approaches the uncorrected wavefront that is dominated by tilt, leading to an 
approximately quadratic surface for the mean square residual error, as seen in Figure 6b. 

As the noise is reduced, the correction becomes more effective, reducing but not 
eliminating the tilt. The shape of the residual error, as shown in Figure 6d, is still 
predominantly quadratic, but the magnitude has been reduced. (Note that the scale is 
varied from case to case in Figure 6 . ) 

In the limit of zero noise, the low spatial frequency terms are eliminated and the 
high-frequency terms dominate in the residual wavefront as in Figure 6h. 
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Correction Function Mean Square Residual Wavefront 




Figure 6. Performance versus noise level. 
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Figure 7 shows the effect on the mean square residual error of decreasing the subaperture 
size. The remaining error is predominantly at the spatial frequency of the sensor array, 
with magnitude proportional to L 5 ' 3 . This fitting error associated with the sensors is 9 
similar to the fitting error for an array of actuators that has been analyzed by Hudgin. 

The overall performance of the first configuration is as a function of noise level and 
subaperture size is shown by the solid curves in Figure 8. For photon noise levels above 
unity, the error is independent of the number of subapertures across the aperture, since 
at these levels only the tilt term is estimated with any reliability. 

At low noise levels, the error depends only on the number of measurements. For a sub- 
aperture dimension equal to the side of the full aperture, the first configuration makes 
two measurements of x and two of y compared to one each for the second configuration. As 
the subaperture dimension is reduced, these edge effects become insignificant and the 
performance of the two configurations becomes identical. 

The apparent difference in the performance of these configurations indicated by the 
analyses of Hudgin 3 and Fried4 has been resolved by including the statistics of the wave- 
front phase as well as the photon noise in the treatment. 

Figure 9 shows the performance for the first and third configurations. When the sub- 
aperture dimension is equal to the side of the aperture, the third configuration consists 
only of the four corners of the array, which is the same as a two-by-two array in the 
second configuration. As in the previous case, the edge and corner effects decrease as 
the subaperture size is reduced. Again, the differences in performance indicated in 
Southwells analysis 8 disappear with a more realistic mpdel and use of a priori knowledge 
of the wavefront statistics. 

Mean Square Residual Wavefront 




Figure 7. Performance versus subaperture size. 
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In adaptive optical systems that compensate for random wave-front disturbances, a wave front is measured and 
corrections are made to bring it to the desired shape. For most systems of this type, the local wave-front slope is 
first measured, the wave front is next reconstructed from the slope, and a correction is then fitted to the recon- 
structed wave front. Here a more realistic model of the wave-front measurements is used than in the previous liter- 
ature, and wave-front estimation and correction are analyzed as a unified process rather than being treated as sepa- 
rate and independent processes. The optimum control law is derived for an arbitrary array of slope sensors and 
an arbitrary array of correctors. Application of this law is shown to produce improved results with noisy measure- 
ments. The residual error is shown to depend directly on the density of the slope measurements, but the sensitivity 
to the precise location of the measurements that was indicated in the earlier literature is not observed. 



INTRODUCTION 

The wave-front correcting system treated here consists of a 
number of sensors, each of which measures the weighted av- 
erage of wave-front slope over some region; a number of cor- 
rector actuators, each of which adds a corrective function to 
the wave-front phase; and a control law relating actuator 
displacements to sensor measurements. The correction will 
be imperfect because of the finite size and number of sensors, 
measurement noise (including photon noise in the sensors), 
the finite number of degrees of freedom in the corrector 
(whether distributed by zones or by modes), and possible 
deficiencies in the control law. 

The approach taken to the control law in the literature has 
been first to reconstruct an estimate of the wave front from 
the slope data and then to fit the actuator functions to the 
estimate. 1 A number of papers have dealt with the problem 
of wave-front reconstruction from noisy data. 2 - 10 These 
papers generally treated the wave front as an array of discrete 
points and the measurements as phase differences between 
pairs of points, and only Ref. 10 included the effect of inte- 
grating wave-front slope over a finite area. The least-mean- 
squares estimate of the wave-front phase based on noisy 
measurements was derived, but no use was made in the pro- 
cess of any knowledge of the statistics of the wave front being 
measured. The latter information is useful even in the case 
of small measurement errors since there are always wave-front 
modes to which the sensor array is blind. 

The problem of fitting a corrector to a wave front has also 
been analyzed by several authors, with particular application 
to wave fronts with the Kolmogorov spectrum typical of at- 
mospheric turbulence. 11 " 15 

In practice, all these problems are interdependent, and the 
derivation of the optimal control law requires a unified 
treatment of the overall system. The present paper gives such 
a treatment for an arbitrary array of wave- front sensors, an 
arbitrary array of corrector functions, and an arbitrary dis- 
tribution of uncorrected wave fronts. The control law that 
minimizes the mean-square residual error and the error itself 
are evaluated. 



SYSTEM DESCRIPTION AND ASSUMPTIONS 

The overall correction system may be represented as a number 
of sensors that measure wave-front slope in the plane of an 
aperture, a corrector with a number of actuators acting in the 
same plane, and a control law connecting sensors to actuators. 
The aperture of the optical system can be described by a 
weighting function proportional to the intensity of light at the 
wave front. It is convenient to normalize this function so 
that 

£~dxW A (x) = l, (l) 

where x is the two-dimensional vector position in the aperture 
plane, W A (x) is the aperture weighting function, and f dx 
indicates integration over the entire aperture plane. 

The uncorrected wave-front phase or optical path difference 
is described by the function i^(x). Because any uniform phase 
over the aperture does not affect optical performance, we can 
equivalently deal with the phase with aperture average re- 
moved, <£(x). For this function 

£dx^(x)*(x) = 0, (2) 

where 

*(x) - rP(x) - dx' W A (x') 

and ^(x) is the input wave-front phase. 

The output of each wave -front sensor is taken as the 
weighted average of the wave-front slope over some region of 
the aperture, which approximates the output of real wave- 
front slope sensors for small wave-front deviations and a fixed 
wave-front intensity pattern. The noise on the measurements 
is represented by a random additive function of position, 
which permits a convenient treatment of photon noise. The 
signal out of the nth sensor is then 

Sn = £~ dx W sn (x)[<t> s n (x) + u(x)] } (3) 



413 



where 

s n is the signal from the nth sensor, 
W sn (x) is the weighting function for the nth sensor, 
<p s n (x) is the slope of wave front in the direction of sensitivity 
of nth sensor, and 

v(x) is the noise distribution over the aperture. 

The signal can be expressed in terms of the wave-front phase 
by integrating Eq. (3) by parts: 



(4) 



where W $ sn {x) is the derivative of W sn (x) in the direction of 
slope sensitivity of the nth sensor. In Eq. (4) use has been 
made of the fact that W sn (x) goes to zero at infinity. 

The noise is assumed to be unbiased and uncorrected with 
the wave-front phase: 



(v(x)) = <y(x)<MxO> =0. 



(5) 



The control law for the system generates a command to each 
actuator of the wave-front corrector based on all the sensor 
outputs. For a linear control law, the command to the ;th 
actuator is 



<7 = T,M jn s n , 



(6) 



where cj is the command to the ;th actuator and Mj n is the 
weighting of nth sensor signal in ;th actuator command. 

Finally, the responses of the actuators are assumed to 
combine linearly to form the total wave-front correction, 



*W = Z cjrjlx), 



(7) 



where <£(x) is the total wave-front correction and r ; -(x) is the 
response of the ; th actuator to a unit command. 

EVALUATION OF MEAN-SQUARE RESIDUAL 
ERROR 

The residual wave-front error may be expressed by using Eq. 
(7) and expanded by using Eqs. (4) and (6): 

€(x) = *(z) - <Mx) 



(8) 



The mean-square residual error at x is formed by using Eq. 
(8) and assuming that the sensor noise and the wave-front 
phase are uncorrected: 

<e2(x)> = IZS Zrj{x)rj>(x)Mj n Mj< n >(s n s n <) 

j y n n' 

-2ZZrj(x)M jn (s n <t>(x)) + <* 2 (x)>, (9) 

/ n 

where (/> is the expected value of /. 

The overall measure of performance is the expected 
mean-square error averaged over the aperture: 



(e?)= f_"dxW A (x)<e*(x)> 



s IEI%M ;V (Mn') 
j j' n n' 



X dxW,(x)r ; (x)r ; v(x) 

-2ZZ Mjn f" dxW A ix)rj(x){s n <p(x)) 

j n 

+ f " dxW A (x)(4> 2 (x)). (11) 

%J ~- CD 

These expressions may be simplified by defining additional 
terms. The integral involving the product of two sensor 
functions is defined as the matrix S nn <: 

= £~ dx' J"_" dx»[Wi B (xW: B (**)<*(x')*(*")> 



+ W sn (x')W sn ,(x"){v(.x')v{x>'))}. 



(12) 



The integral involving the product of two actuator response 
functions is defined as the matrix Rjy: 



R jr = f" AxW A (x)rj{±)r r (x). 



(13) 



The integral involving products of sensor and actuator re- 
sponse functions is defined as the matrix Aj n : 

A jn = £° dxW A (x)r,(x)<s„tf>(x)) 

= f°dx f° dx'W /4 (x) O ( X )W s sn (x')(0(x)^(x')>. 

%J —co *J — co 

(14) 

The average mean-square uncorrected error is defined as 

(4): 

= J > _°dx^(x)(02(x)>. (15) 

Substituting these definitions into Eq. (11) gives the simpli- 
fied form 



(4) = E E E E M jn M rn >s n „, Rjr -2EE M jn A jn + (4). 

j j' n n' j n 



(10) 



(16) 



MINIMIZATION OF RESIDUAL ERROR 

Equation (16) gives the average mean-square residual error 
for arbitrary weights in the control matrix Mj n . To optimize 
system performance we wish to choose the weights that min- 
imize the error. 

Differentiating Eq. (16) with respect to element My n > and 
equating to zero yields 

7T7^ " 2 £ £ M jn S nn <Rjj< - 2A rn > 

OMj'n' j n 

= 0. (17) 

By letting repeated indices indicate summation and making 
use of the symmetry of %<, Eq. (17) may be written more 
concisely as 

Rj'jMjnSnn' = A/ V . (18) 
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The solution to this equation is the minimizing control matrix 

M'jn'' 

M) n = RpAjwS^ (19) 

The R matrix will not be singular unless there are redundant 
actuators that allow the same correction to be made with 
different combinations of actuator commands. In that case 
the generalized inverse of the matrix may be used at the ex- 
pense of additional computation, or the redundant functions 
can be eliminated with equivalent results. 

Alternatively, a cost function depending on the actuator 
commands can be added to the quantity being minimized. If 
the cost rises more than linearly with the actuator command, 
the correction will be distributed over the redundant modes, 
and R will not be singular. The S matrix will not be singular, 
even with redundant sensors, if there is independent noise on 
each measurement, as will be the case in any real system. 

The value of the minimum average mean-square residual 
error, <€?>, is found by substituting Eq. (19) in Eq. (16): 



(£) = (el) - Mj n A jn 

= (4) - RtfAjvSSAjn. 



(20) 



This expression, with Eqs. (12)-(14), gives the optimal per- 
formance for an arbitrary system and arbitrary wave-front 
statistics. 

The overall performance will depend on the magnitude of 
the sensor noise relative to the expected signal that is due to 
fluctuations of the input wave front. If the noises in all sen- 
sors, <rj, are equal and are large compared with the signal 
variance, the matrix approaches //<r§, the correction term 
in Eq. (20) becomes small, and the residual error approaches 
the uncorrected error. 

For small measurement errors, the residual error will ap- 
proach a limit set by the inability of the finite arrays of sensors 
and actuators to fit the wave front perfectly. 

In the transition region where the residual error is small 
compared with the uncorrected error but large compared with 
the fitting error, the performance will depend directly on the 
measurement errors and the error propagation of the control 
matrix. 



The correlation required is 

<0(x)0(xO> - dx" J_~ dx"'W A (x")W A (x"') 

= -V 2 f~dx" C dx"'W A {x")W A {x"') 

%} CO %J —CO 

X [D(x, x') - £(x, x'") - D(x", x') 
+ D(x",x")] 

= -V 2 D(x, x') + g(x) + *(x') - a, (23) 

where 

£(x) = V 2 f'dx'W^x'JDfc.x*), (24) 

•/ —CO 

a = V 2 §~ dx" J~ dx"W A (x")W A (x")D(x",xn. 

(25) 

Substituting into Eq. (15), 

Ul) = J^dxW A (x)(4>Hx)) 

- a. (26) 

In the evaluation of S and A, all terms on the right-hand side 
of Eq. (23) except the first can be seen to integrate to zero, 
leaving 

Snn>= C dx' f ~ dx*[-y 2 w\ n (x')w\A*")D(x\xn 

%J — OO xJ — CO 

+ W sn (x')W sn <(x"Hv(x')v(x"))], (27) 
Aj n = -y 2 £~ dx J~ dx'W A (x) rj (x)W s sn (x')D(x, x'). 

(28) 

Equations (26)-(28) may be used with Eqs. (19) and (20) to 
express the optimal control law and evaluate its performance 
in terms of the structure function of the wave front being 
corrected. 



REPRESENTATION BY STRUCTURE 
FUNCTION 

In the case of atmospheric turbulence, the statistics of the 
input wave front are defined in terms of a structure func- 
tion, 



D(x,x0= -*(x')F>. 



(21) 



The correlation of phases with averages removed required in 
the computation of S, A, and (el) can be derived directly by 
using Eq. (21). The average removed phase is 

<i>(x) = fix) - f~ dx"W*(x")lMx") 

= f^dx"W A (x"M(x) (22) 



PHOTON NOISE 

In many systems of interest, the limit on performance is set 
by the noise that is due to photon statistics in the wave-front 
measurements. The photon noise for a single wave front can 
be represented by making u(x) a spatially white-noise func- 
tion. The correlation of the noise on two wave-front mea- 
surements depends in detail on the configuration of the 
wave- front sensor used. If the measurements are made on 
separate wave fronts, either in sequence or separated by a 
beam splitter, the noises will be uncorrected, even though 
they refer to the same area in the aperture. 

Measurements of orthogonal slopes made simultaneously 
on the same wave front will also have uncorrelated noises. To 
reflect these conditions the correlation of the noise functions 
may be written as 



(v(x)v(x')) =/w<r 2 n 5(x-x'), 



(29) 
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where 

<r 2 is the photon noise density, 

k n n' m In * V for n and n' measured simultaneously on the 
same area of the same wave front but 
knn' is 0 otherwise, 

!„ is a unit vector in the direction of slope sensitivity of the 
nth measurement, and 

5(x) is the Dirac delta function. 

The noise term in Eq. 10 or 27 may then be written as 

f " dx' f " dx'WsnWWsnix'Hvix'M*")) 

%J — CO — CD 



= k ,<T 2 



W sn (x')W m >(x'). 



(30) 



In most cases the basic measurements will be disjoint, and 
h nn > will be the Kronecker delta if each measurement is treated 
separately. It may be convenient to combine a number of the 
basic measurements before processing in order to reduce the 
order of the matrices involved in the computation of the 
control law. In this case the more general form of k nn ' would 
be needed. 

CLOSED-LOOP RECONSTRUCTION 

Most wave-front correction systems that have been imple- 
mented up to the present time have used a feedback-control 
loop in which the actuators are driven to null the wave-front 
sensor signals generated by the input wave front. Such a 
system is analyzed here for comparison with the optimum 
system derived above. 

If a set of commands cy is sent to the actuators, the phase 
correction created will be given by Eq. (7). Letting repeated 
indices indicate summation, we may write this as 



4>(x) = cjrj(x). 



(31) 



If this correction is subtracted from the input wave front, the 
signal from the nth sensor as given in Eq. (4) becomes 

s n '=s n - f_^dx[-Wl n (x)4>(x)] 



(32) 



The integral on the right-hand side of Eq. (32) for a unit ac- 
tuator command defines the matrix P n j, which relates the nth 
sensor output to the jth actuator input: 



Pnj « £ dx[-W sn (x)rj(x)}. 



(33) 



Substituting in Eq. (32), we find that 

s n =s n - P nJ Cj. (34) 

In the closed-loop operation considered here, the command 
vector cj is chosen to drive the corrected signal vector s' n to 
zero. This requires that 

Cj=PijS n , (35) 

where PJ) is the generalized inverse of P n j . The generalized 
inverse minimizes the sum of the squares of the residual sig- 



nals and sets to zero any actuator modes not detected by the 
sensors. 

The expected mean-square error averaged over the aperture 
is found by substituting the control matrix Pj n for M Jn in Eq. 
(16): 



(6?) = (el) - 2P + nj A Jn + PZjP + nr S nn >Rjj, 



(36) 



For a system with a given configuration of sensors and actu- 
ators, Eq. (36) allows the performance using simple closed- 
loop control to be compared with that obtained using optimal 
control. 



COMPUTATIONAL RESULTS 

The theory derived here has been applied to a system con- 
sisting of a square aperture with uniform illumination and a 
square array of identical actuators with Gaussian response 
functions. The array is positioned with an actuator at each 
corner of the aperture, and the ratio of the side of the aperture 
to the actuator spacing, N, is varied in the computations. The 
actuator response falls to lie at the adjacent actuator. 

The wave front being corrected is assumed to have the 
Kolmogorov structure function typifying atmospheric tur- 
bulence: 

D(x, x') = 6.8839pr - ^O/r^/a, (37 ) 

where D{x, x') is the phase-structure function, rad 2 , and r 0 
is Fried's coherence length, m. The mean-square uncorrected 
phase error over the full aperture for this structure function 
is 



(e 2 Q ) = 1.3103(A/r 0 ) 5/3 , 



(38) 



where (el) is the uncorrected phase error, rad 2 , and A is aside 
of aperture m. 

The wave-front sensor investigated averages the wave-front 
slope with uniform weight over a square subaperture of di- 
mensions equal to the actuator spacing. Two different con- 
figurations of such sensors are considered. 

In the first configuration, shown in Fig. 1, X and Y slopes 
are measured in displaced subapertures. This configuration 
is typically used in the shearing interferometer sensor and was 
treated by Rimmer 2 and Hudgin. 4 The mean-square slope 
measured over a square subaperture of side L for the Kol- 
mogorov structure function is 



c] = 6.40511-^-5/^ 



(39) 



where <r 2 is the mean-square slope signal, rad 2 m -2 , and L = 
A/N is a subaperture dimension m. The noise on this mea- 
surement is assumed to be pure photon noise for which the 



Aperture 



Actuator 

location 



Sensor 
Subaperture 



Fig. 1. Displaced subaperture wave-front sensor configuration. 
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mean-square slope error is inversely proportional- to the 
number of photons used in the measurements. 

In this configuration, a sensor subaperture located at an 
edge of the aperture parallel to the slope being measured will 
have only half of its area illuminated and will therefore have 
twice the mean-square slope error of a fully illuminated 
subaperture. 

For a uniform illumination yielding p photons per square 
meter, the noise will be 

<Tn M pA s ' 

where 

a* is the mean-square ,slope noise, rad 2 m~ 2 , 
K is a constant, rad 2 m~ 2 , 
At is the number of photoelectrons used, 
p is the photoelectron density m~ 2 , and 
A s is the subaperture area m 2 . 

Figure 2 shows the residual error for this configuration as a 
function of photoelectron density and actuator spacing. The 
aperture size A and r 0 , where each is assumed to be 1 m, and 
a photon noise constant of 1204.6 were used. (This constant 
would apply to a shearing interferometer with optimum shear 
imaging a disk of 5-^rad diameter in light of 0.55-^m wave- 
length.) 

The three regions of operation are easily distinguished in 
-Fig. 2. For densities much less than 100, the measurement 
noise is so large that the a priori estimate of zero phase is more 
accurate than the measurement For high densities, the error 
is determined by actuator and sensor spacing and is again 
independent of irradiance. This fitting error varies as (A/ 
Nro) b/s when N is sufficiently large that edge effects may be 
neglected. 

In the transition region, where the photon error is small 
compared with the uncorrected error but large compared with 
the fitting error, the performance depends primarily on 
photoelectron density and depends relatively weakly on ac- 
tuator spacing, in agreement with the results of Fried 3 and 
Hudgin. 4 

Figure 3 compares the performance of the optimal recon- 
structor with that of the closed -loop reconstructor for N - 4. 
At low photoelectron densities, the closed-loop reconstructor 
faithfully follows the noise input, leading to large errors. At 
high densities, the performance is close to that of the opti- 
mum, but the actuators are not fitted to the wave-front quite 
as well as possible. The optimal reconstructor thus improves 
performance at all illumination levels. 

The second sensor configuration considered is shown in 
Fig. 4. Here the x and y slopes are measured in a common 
subaperture that has an actuator at each corner. This con- 
figuration characterizes the Hartmann sensor and was ana- 
lyzed by Fried. 3 

The performance of the two configurations is compared in 
Fig. 5. The top dashed curve is for a single subaperture with 
four actuators, one at each corner. The fitting error in this 
case is slightly larger than for the displaced subaperture 
configuration with the same subaperture size, which uses twice 
as many slope measurements to link the four actuators. For 
larger numbers of actuators for which edge effects are reduced, 
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Fig. 2. Performance of displaced subaperture sensor configura- 
tion. 
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10 
5 

2 

t 1 
J °* 5 
f 0.2 

! ot 

i 

0.05 
0.02 









A * 1 m 

r Q * 1 m 

a 2 • l,204.6/( 




















1 1 — 


- Displaced suba 

- Combo n suba per 

, ■ 


aerture eonfigur 
ture configurati 


ation 
on 





10 1 10 : 10 5 10- 10* 10* 

Fig. 5. Comparison of sensor configurations. 

the performance of the two configurations is virtually indis- 
tinguishable. The difference in performance between these 
configurations as given by Fried 3 and Hudgin 4 does not appear 
in these results. 
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As was discussed by Herrmann, 8 when the wave-front re- 
construction algorithm based on the average of the phases at 
the four nearest neighbors plus the average of the four con- 
necting phase differences is used with the common subap- 
erture configuration, two interleaved arrays of points are 
produced with an arbitrary phase difference between them. 
Fundamentally, the approach given here resolves the ambi- 
guity by using the statistics of the wave fronts being estimated 
to imply the unobserved mode. 

These results also contradict those of Southwell, 9 who shows 
a mean-square error for the common subaperture configura- 
tion approximately twice that of the displaced subaperture 
configuration. In his model, slopes were measured at discrete 
points with equal noise on all measurements, and the wave- 
front error was evaluated only over an array of discrete 
points. 

With more-realistic sensors that average the slopes over 
finite subapertures, with a performance measure that evalu- 
ates the error over all the aperture, and with use of the sta- 
tistics of the wave-front phase in the estimation process, the 
extreme sensitivity of performance to the precise sensor 
configuration reported by Southwell is not observed. 

The same conclusion was reached in a separate paper 16 that 
dealt with the same sensor configurations but allowed actuator 
response functions to be chosen arbitrarily. Again, except for 
edge effects, performance was the same for the two configu- 
rations treated here as well as the configuration considered 
by Southwell in which X and Y slopes are measured in com- 
mon subapertures centered at the actuator locations. 

The results of the performance computations indicate that 
appreciable improvement in performance can be obtained by 
using the optimal reconstructor and that performance de- 
pends on sensor dimensions but not significantly on the po- 
sition of the X or Y sensor array with respect to the actuator 
array. 
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Lawrence W* Stark dies at 78 





Related links: 

UC_^ress_reiease 
leieroboticsJuTeiiroi^ 
Conversations with history 



Dr. Lawrence W. Stark, a professor emeritus of physiological 
optics and engineering at the University of California, Berkeley, 
recognized worldwide as a pioneer in the use of control and 
information theory to characterize neurological systems, died 
Friday, Oct. 22. 



Stark died of cancer at his home in Berkeley at the age of 78, 
four years after being diagnosed with Non-Hodgkin's lymphoma. 



"Stark's death is a huge loss for the world of optometry," said 
Dennis Levi, dean of UC Berkeley's School of Optometry. "His 
landmark studies on eye movement control truly advanced the 
field of vision science." 



A neurologist by training, Stark is credited for seminal research 
that applied engineering principles, particularly control theory, 
to biological systems. 



Return to Anmuricements 



"Stark was unique in his ability to identify aspects of 
engineering analysis relevant to medicine and biology," said 
Laurence Young, Apollo Program professor of astronautics at 
the Massachusetts Institute of Technology and Stark's first 
graduate student at MIT. "A perfect example is the seminal 
work he conducted in analyzing the way a pupil reacts to light in 
terms of a linear control system." 



Young pointed out that those same principles have been applied 
to such areas as pilot control of airplanes and spacecraft. 
"Characterizing how the eye moves and how the brain 
processes visual cues is essential to understanding how pilots 
control airplanes, and why people get motion sickness," he said. 



In addition to analyzing the feedback control system governing 
pupil contractions, Stark also developed the scan path theory of 
eye movements. He studied the way people's brains viewed the 
world and analyzed the vast number of jumping eye 
movements, or "saccades," people make. He noticed specific 
sequences to how people glimpsed a room, face or other scene 
before them, and realized how those sequences provided clues 
to the importance of pictures generated by the brain. 



://www.me.berkeley.edu/announcements/stark.html 
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Stark's drive to understand visual processes within an 
engineering discipline led to his later research interests in 
robotic vision and virtual reality. 



"Although he was trained in medicine, he was very interested in 
the physical sciences and engineering/' said Gerald 
Westheimer, UC Berkeley professor of neurobiology. "He was a 
true crossover scientist, bringing applied engineering concepts 
to neurological functions, and the variables inherent in biology 
to engineering." 



Stark was born in New York on Feb. 21, 1926. He credited his 
early interest in engineering to the influence of his father, an 
MIT-educated chemical engineer. He once recalled in an 
interview how, as a young boy, he took apart his mother's 
typewriter and put it back together again — minus four screws. 
His mother was impressed with his success, he said, until the 
typewriter stopped working a few weeks later. 



Undaunted, Stark maintained his curiosity for how things work, 
going on to Columbia University and majoring in English, 
biology and zoology. After receiving his bachelor's degree in 
1945, he joined the military, taking the U.S. Navy up on its 
offer to pay for his medical school tuition. Stark went to New 
York's Albany Medical College, where he earned his M.D. in 
1948. 



He then spent two years in England at Oxford University and 
University College, conducting research in neurophysiology, 
biochemistry and biophysics, before returning to the U.S. Navy 
to serve as a doctor during the Korean War. 



In 1954, after the war ended, Stark joined Yale University as an 
assistant professor of medicine. 



In 1960, he became head of the neurology section of MIT's 
Center for Communication Sciences, and in 1965, he founded 
and became chairman of the Biomedical Engineering 
Department, one of the country's first bioengineering 
departments, at the University of Illinois at Chicago. 



"He was really one of the first people ever to use engineering 
theory to study a physiological system," said Blake Hannaford, 
director of the Biorobotics Laboratory at the University of 
Washington and one of Stark's former Ph.D. students. "He 
played a pivotal role in the 1950s and 1960s in founding the 
field of bioengineering." 
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Stark joined the UC Berkeley faculty in 1968 with joint 
appointments as professor at the School of Optometry, the 
Department of Electrical Engineering and Computer Sciences 
and the Department of Mechanical Engineering. He also 
collaborated on neuro-ophthalmology research with colleagues 
at UC San Francisco. 



He retired from UC Berkeley in 1994, but remained active in his 
lab on campus. 



"He was always questioning, bringing up new ideas/' said Elwin 
Marg, UC Berkeley professor emeritus of vision science and 
optometry and a close colleague of Stark's for nearly 50 years. 
"His enthusiasm and curiosity were an inspiration to his 
students." 



So many of his students went on to distinguish themselves in 
academic careers that a scientific conference was held in his 
honor. In 1994, John Semmlow, a professor of bioengineering 
at Rutgers University in New Jersey and one of Stark's former 
students, organized the First International Conference on Vision 
and Movement in Man and Machine, a two-day symposium 
attended primarily by Stark's colleagues and former students 
and affectionately nicknamed "Starkfest." 



Although research papers were presented at the conference and 
later published in peer-reviewed journals, including special 
issues of the Annals of Biomedical Engineering and Optometry 
and Vision Science, participants viewed the event as a chance 
to honor their former advisor and renew old friendships, said 
Semmlow. 



"Stark maintained a strong bond with nearly all of his students," 
said Semmlow, one of two students from the University of 
Illinois who followed Stark to UC Berkeley. "He had an 
extremely dynamic personality. He was also extraordinarily 
intelligent, very well read, and he cast his interests in many 
different directions." 



Two more Starkfest meetings have been held since 1994, most 
recently in Marseilles, France, in 2002. 



Stark received numerous honors throughout his career, 
including an honorary Sc.D. from the State University of New 
York and an honorary Ph.D. from Tokushima University in 
Japan. He was named a Guggenheim Fellow in 1968, a fellow of 
the Institute of Electrical and Electronics Engineers in 1970, a 
recipient of the William J. Morlock Award in Biomedical 
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Engineering in 1977, and a fellow of the American Institute of 
Medical and Biological Engineering in 1992. 



Stark is survived by his partner of 18 years, Jill Strohn of 
Berkeley, Calif.; three daughters, Stefanie Stark of Kensington, 
Calif., Nanou Matteson of Berkeley, Calif., and Elizabeth Stark 
of San Francisco; Elizabeth's mother, Wendy Bartlett of 
Berkeley, Calif.; ex-wife, Jeanne Stark-Iochmans of Berkeley, 
Calif.; his brother, Matthew of Minneapolis, Minn.; and four 
grandchildren. 



Stark's first wife and Stefanie's mother, Jane Stark, died in 
2001. 



Memorial services at the UC Berkeley Faculty Club are being 
planned. When the exact date and time are available, they will 
be posted online at http://scan.berkeley.edu/larry . In lieu of 
flowers, donations in Stark's memory can be made to the 
Golden Gate National Parks Conservancy, Fund in Memory of 
Lawrence W. Stark, Attn: Audrey Yee, Fort Mason Building 201, 
San Francisco, CA 94123-0022. In accordance with Stark's 
wishes, the money will be used to purchase land that will be 
kept preserved and open to the public 
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